Numerical integration - Trapz, quad, quadl...
in [Matlab]
|
From: DD C on 9 Feb 2010 11:22 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!
From: Cygnine on 9 Feb 2010 11:31 "DD C" <dchr0036(a)remove.thisgmail.com> wrote in message > 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.
From: DD C on 9 Feb 2010 14:10 "Cygnine " <cygnine(a)remove.this.gmail.com> wrote in message <hks2k9$bts$1(a)fred.mathworks.com>... > 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: ============== >> 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 ============== >> F = @(x)(1/0.0079)*exp(-50000.5*x).*x.^(0.5).*(1-x).^(64.5); >> 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.
From: Cygnine on 9 Feb 2010 14:55 "DD C" <dchr0036(a)remove.thisgmail.com> wrote in message > 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 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.
From: Walter Roberson on 9 Feb 2010 15:03
DD C wrote: > 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: > ============== >>> 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 |