Numerical integration - Trapz, quad, quadl...

[DD C]

Hi,

I am trying to calculate several integrals. It appears that in some cases, trapz provides a result while quad and quadl provides 0. In other cases quad and quadl provides results while trapz returns NaN.
I understand the approximations behind each method but I don't know why one function provides results while the other doesn't. Are there specific cases when quad, quadl, trapz should not be used? or Are there specific cases when they should be used?

Thanks!

[Cygnine]

Post by DD C
I am trying to calculate several integrals. It appears that in some cases, trapz provides a result while quad and quadl provides 0. In other cases quad and quadl provides results while trapz returns NaN.

Can you provide an example that reproduces this behavior? Trapz works on an array, right? But quad and quadl work on functions...so how do you give data to trapz? With equispaced function evaluations? What is the function you're trying to integrate?

Off the top of my head, I don't see why quad (or quadl) would be worse than trapz unless you're trying to integrate something that's quite pathological or not integrable.

[DD C]

Post by Cygnine
Can you provide an example that reproduces this behavior? Trapz works on an array, right? But quad and quadl work on functions...so how do you give data to trapz? With equispaced function evaluations? What is the function you're trying to integrate?
Off the top of my head, I don't see why quad (or quadl) would be worse than trapz unless you're trying to integrate something that's quite pathological or not integrable.

I am trying to integrate the following function between 0 and 1:
f = (1/0.0079)*exp(-50000.5000*x)*x^(0.5)*(1-x)^64.5

I tried the following two approaches:
==============

Post by Cygnine

x = 0:0.00001:1;
y = (1/0.0079)*exp(-50000.5*x).*x.^(0.5).*(1-x).^(64.5);
z = trapz(x,y)

ANS
z = 9.2370e-06
==============

Post by Cygnine

Q = quad(F,0,1)

ANS
Q = 0
==============

I suspect it comes from the fact that the answer is a very small number. I also tried to play around with the tolerance value of the "quad" function but I am not able to produce another number than 0 with "quad".

Thank you for your help.

[Cygnine]

Post by DD C
f = (1/0.0079)*exp(-50000.5000*x)*x^(0.5)*(1-x)^64.5
I suspect it comes from the fact that the answer is a very small number.

That's...quite a function you have there. First, if you use double precision arithmetic (default in Matlab), you'll never compute anything other than zero for the integral over [1.5e-2, 1]. The function value is below the minimum absolute value that's resolvable (type "realmin").

So you cannot compute an integral directly over [1.5e-2, 1] using any direct quadrature method. (You may be able to do some funky transformations with logs if you play around with the math for a bit.)
However, if you try
quad(F, 0, 1.5e-2)
you'll see that you do get some nonzero answer. Which one is more accurate (quad, quadl, trapz)...who knows; that's a pretty crazy function you're trying to integrate. I would break it up into small subintervals (width 1e-5 or so) and iterate up until 1e-2 or whatever. For example,
quad(F, 0, 1e-5)
and
quadl(F,0, 1e-5)
agree to a few digits.

[Walter Roberson]

Post by DD C
f = (1/0.0079)*exp(-50000.5000*x)*x^(0.5)*(1-x)^64.5
==============

x = 0:0.00001:1;
y = (1/0.0079)*exp(-50000.5*x).*x.^(0.5).*(1-x).^(64.5);
z = trapz(x,y)

ANS
z = 9.2370e-06

FYI, the symbolic solution to a number of decimal places is 0.00001001420547.
The accuracy of this will be much affected by the accuracy of the 1/0.0079
If we replace the 1/0.0079 with 1/p and solve symbolically, we get

5/12426501021958097388107759726327813799660233187055127117518234795912902964900750923555951686456645680527245612555012624593718416824382754263347392890529665812406156646440027875303740166659395396137811962369232259
8653321549212879464735126538420000481410506974156454636951175709562236623419147958877609548578322146659132902
* Pi * exp(-100001/4) *
(124025561061247309457529144460208287306426791146246195331111079177506959961017477858779094376652204917243208996349070422256031615300455041283710367073742211568707518015092680651202429554214702\
263900182631704459898851239029090474931300725941518249117237902866526465479470169183278534599947731784542673739499452723055
* BesselI(0,100001/4) +
124023080649241820659374163695527605198214470628588899709418432371437299966478305321588018683485392915198176681754682611995195541468339367737626270850133274380181533563254951330\
905648302553861938551651571095765736369724189459763237709303350504853292346637877691117509551064142702152568167094217355641718433763059427
* BesselI(1,100001/4)) / p

[John D'Errico]

Post by DD C
Hi,
I am trying to calculate several integrals. It appears that in some cases, trapz provides a result while quad and quadl provides 0. In other cases quad and quadl provides results while trapz returns NaN.
I understand the approximations behind each method but I don't know why one function provides results while the other doesn't. Are there specific cases when quad, quadl, trapz should not be used? or Are there specific cases when they should be used?

I might argue that if you truly do understand the
methods involved, then you will also understand their
limitations, as well as guidelines where you might use
one over another.

Trapz is a simple tool that integrates a function
represented as a fixed sequence of function
evaluations. So you already have a data series here,
and you need to know the integral of the function
this series is sampled from.

Trapz might be applied to data that you have
measured, so there might be noise in these
measurements. Note that when noise is present,
a lower order method (such as trapz) is a good
choice, as it will minimize the variance of the
estimated integral due to that noise.

The quad family of tools are ADAPTIVE methods.
They are useful when you have a complete black box
of a function to integrate. These tools require the
ability to evaluate the function at a programmatically
specified set of points, where matlab chooses the
specific sequence for these evaluations.

Quad and its cousins can be problematic when you
have a function with a spike in it. (A spike is a
delta function, or something near to that shape.)

As such, the quad tools are very distinct from trapz.
You will want to use them at different times, based
on your goals, on what form your "Function" is in, etc.

John