Prev: How to fit with an infinite serie function ?
Next: Embedded MATLAB Function
From: leo nidas on 9 Jul 2010 06:09

Hi there,

I am trying to use quadgk to numerically evaluate an integral. I have shown that the following integtral (K) is infimum for the values I have selected in the second line. Indeed if I use "int" I take "Inf" as an answer. For the same intergal if I use quadgk I use 0.0580.

The worrying fact is that quadgk yields the result without a warning. (Matlab 2007b).
Why is that?

syms x real
a=0.8;t=0.2;g=2;l=3;
fx=exp(x.*t-g.*l.^a.*x.^a).*x.^(a.*g-1);
K=int(fx,'x',0,Inf)

fx=@(x) exp(x.*t-g.*l.^a.*x.^a).*x.^(a.*g-1);
K=quadgk(fx,0,Inf)
From: Steven Lord on 9 Jul 2010 10:39

"leo nidas" <bleonidas25(a)yahoo.gr> wrote in message
news:i16sg2$j07$1(a)fred.mathworks.com...
>
> Hi there,
>
> I am trying to use quadgk to numerically evaluate an integral. I have
> shown that the following integtral (K) is infimum for the values I have
> selected in the second line. Indeed if I use "int" I take "Inf" as an
> answer. For the same intergal if I use quadgk I use 0.0580.
> The worrying fact is that quadgk yields the result without a warning.
> (Matlab 2007b).
> Why is that?

For small values of x, the value of fx is very small. But somewhere between
10^6 and 10^7, the value of the integrand shoots up to become very large.
The symbolic integration routine probably checks the limit of the integrand
as x approaches Inf (to make sure it's actually approaching 0) and because
it does not returns Inf. It can do this because it can calculate the limit.

But the numeric integration routine performs calculations in double
precision and it needs to choose a specific set of points at which to
evaluate your function. If you write a function that displays the points
QUADGK chose and the function value (IN DOUBLE PRECISION) of your fx at
those points, you'll find that almost all of those values are very, very
small. If x is small then the value of the exponential is small and so is
the value of x^(3/5), the term by which the exponential is multiplied.

If x is large (but not greater than the threshold where the integrand
skyrockets) then the value of the exponential underflows to exactly 0
because the expression inside the exp() call becomes a large magnitude
negative number. Thus the value of the whole expression is 0. QUADGK
evaluates the function at enough points to make QUADGK think that the
function value past a certain point is a constant 0 and so it doesn't need
to refine its estimate, and then it only focuses on the regions where it
received a positive value for the function.

It turns out that by default QUADGK does not choose any of the points at
which it evaluates your function with a large enough value to reach the
point where the value of the integrand skyrockets; if you were to force it
to do so by setting Waypoints, it will warn you that the value of the
function (IN DOUBLE PRECISION) is Inf and QUADGK will return Inf.

The moral of the story is that you shouldn't blindly throw functions and
limits into any of the numeric integration routines -- know your functions
and know where they exhibit "interesting" behavior and tell the integration
routine to focus on those "interesting" regions.

--
Steve Lord
slord(a)mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com


 | 
Pages: 1
Prev: How to fit with an infinite serie function ?
Next: Embedded MATLAB Function