Re: Re: Re: Integrate vs NIntegrate

*To*: mathgroup at smc.vnet.net*Subject*: [mg89211] Re: [mg89182] Re: [mg89147] Re: [mg89129] Integrate vs NIntegrate*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Thu, 29 May 2008 07:05:33 -0400 (EDT)*References*: <9489155.1211802216995.JavaMail.root@m08> <op.ubr1yjh52c6ksp@bobbys-imac> <2b8d8f40805261936u2ec81454r1b6ba16be0897662@mail.gmail.com> <200805271114.HAA01380@smc.vnet.net> <29744585.1211975587467.JavaMail.root@m08> <op.ubvkuekr2c6ksp@bobbys-imac>

Perhaps you should pay a little more attention to other posts in the same thread (not just mine, however flatterning I find this). In case you have not noticed, I strongly recommend David Cantrell's posts, particularly on such kind of topics and, of course, Daniel Lichtblau. Although I really do not see why you need my help in finding them, here are the relevant quotes. Here is David: > Res1 will not be the desired value. Reason: f is not an > antiderivative of h > on the interval [0, 3] due to the fact that f has a discontinuity at > x = 2. > > ........... > > We can, nonetheless, use f to get a reasonably accurate value for the > integral: > > In[20]:= N[f[b] - f[2 + 10^-7] + f[2 - 10^-7] - f[a]] > Out[20]= 1.65267 + 0. I and here is Daniel: > The antiderivative, f[x], has a jump discontinuity due to a branch cut > that is either at or very near to x=2 (probably at). Have a look at > > Plot[{Re[f[x]], Im[f[x]]}, {x, 0, 3}] > > Thus evaluation at the endpoints will fail to recover the definite > integral. Integrate[h[x], {x, a, b}] will give a more plausible > result, in > accord with the NIntegrate value. By the way, I would have been very happy to have staked some real money on this but I feel it would have been, kind of, unfair. Andrzej Kozlowski On 29 May 2008, at 02:43, DrMajorBob wrote: >> Not a proof, of course, but I would bet a fair amount on it >> [Integrate] being right. > > I believe someone recently said we should ignore predictions from > those who haven't staked actual money on it. > > But anyway, is Integrate correct in BOTH these calculations? > > q = -1/2; a = 0; b = 3; > h[x_] = (1 + x^3)^q; > f[x_] = Integrate[h[x], x]; > f[b] - f[a] // N > > -2.55387 - 2.42865 I > > f[x_] = Integrate[h[t], {t, 0, x}]; > f[b] - f[a] // N > > 1.65267 > > Along with the OP, I'm betting the first result from Integrate is > wrong, in view of > > NIntegrate[h[x], {x, a, b}] > > 1.65267 > > If NIntegrate is wrong instead, that hardly improves matters. > > Bobby > > On Wed, 28 May 2008 03:46:27 -0500, Andrzej Kozlowski <akoz at mimuw.edu.pl > > wrote: > >> I would say that Integrate is correct. Evidence: >> >> In[1]:= q = -2^(-1); >> >> In[2]:= f = Integrate[(1 + x^3)^q, x]; >> >> In[3]:= g = Chop[N[D[Normal[f + O[x]^12], x]]] >> >> Out[3]= -0.3125000000000013*x^9 + 0.37500000000000105*x^6 - >> 0.5000000000001258*x^3 + >> 0.9999999999999998 >> >> and, on the other hand: >> >> In[4]:= N[Normal[(1 + x^3)^q + O[x]^10]] >> Out[4]= -0.3125*x^9 + 0.375*x^6 - 0.5*x^3 + 1. >> >> >> Not a proof, of course, but I would bet a fair amount on it being >> right. >> >> Andrzej Kozlowski >> >> >> On 27 May 2008, at 20:14, DrMajorBob wrote: >> >>> Because they're completely different. One is an antiderivative in >>> Complex >>> variables; the other is a definite integral from 0 to x (so you know >>> that >>> it is zero at zero, among other things). Antiderivatives tend to >>> have >>> branch cuts, poles, etc., and Mathematica's chosen branch is complex >>> in >>> part of your region. >>> >>> Or... it's entirely possible Mathematica's answer for >>> Integrate[h[x], x] >>> is simply WRONG. >>> >>> Let's check whether it has the correct derivative: >>> >>> q = -1/2; a = 0; b = 3; >>> h[x_] = (1 + x^3)^q; >>> f[x_] = Integrate[h[x], x]; >>> FullSimplify[f'[x] == h[x], a < x < b] >>> >>> Sqrt[-(1 + x) (-1 + (-1)^(1/3) x)] == ( >>> 2 (-1)^(1/12) Sqrt[-(-1)^(1/6) ((-1)^(2/3) + x) (1 + x^3)])/( >>> Sqrt[1 - I/Sqrt[3] + (-1 - I/Sqrt[3]) x] Sqrt[ >>> 3 + I Sqrt[3] + I (3 I + Sqrt[3]) x]) >>> >>> As you see, FullSimplify doesn't think Integrate got a correct >>> result. Or >>> at least, it couldn't verify it. >>> >>> Here's an attempt to provide FullSimplify some adult supervision: >>> >>> simple= FullSimplify[ #, a<x<b]&; >>> test= f'[x]== h[x]//simple >>> >>> Sqrt[(-(1 + x))*(-1 + (-1)^(1/3)*x)] == >>> (2*(-1)^(1/12)*Sqrt[(-(-1)^(1/6))*((-1)^(2/3) + x)*(1 + x^3)])/ >>> (Sqrt[1 - >>> I/Sqrt[3] + (-1 - I/Sqrt[3])*x]*Sqrt[3 + I*Sqrt[3] + I*(3*I + >>> Sqrt[3])*x]) >>> >>> test2 = simple[test /. {(-1)^(1/3) -> -1, (-1)^(2/3) -> >>> ComplexExpand[Exp[I*(Pi/3)]], (-1)^(1/6) -> ComplexExpand[Exp[I*(Pi/ >>> 6)]], >>> (-1)^(1/12) -> ComplexExpand[Exp[I*(Pi/12)]]}] >>> >>> 1 + x == (((1/2 + I/2)*(3 - I*Sqrt[3])*Sqrt[(-(2*I + (I + >>> Sqrt[3])*x))*(1 >>> + x^3)])/Sqrt[3 - I*Sqrt[3] + (-3 - I*Sqrt[3])*x])*Sqrt[3 + >>> I*Sqrt[3] >>> + I*(3*I + Sqrt[3])*x] >>> >>> test3 = test2 /. (z_)^(r_) :> ComplexExpand[z]^r >>> >>> 1 + x == (((1/2 + I/2)*(3 - I*Sqrt[3])*Sqrt[(-Sqrt[3])*x - >>> Sqrt[3]*x^4 >>> + I*(-2 - x - 2*x^3 - x^4)])/Sqrt[3 - 3*x + I*(-Sqrt[3] - >>> Sqrt[3]*x)])*Sqrt[3 - 3*x + I*(Sqrt[3] + Sqrt[3]*x)] >>> >>> simple[(#1^2 & ) /@ test3] >>> >>> 2*I + (3*I + Sqrt[3])*x == 0 >>> >>> I don't think Integrate has the right answer. >>> >>> Bobby >>> >>> On Mon, 26 May 2008 21:36:31 -0500, Armen Kocharyan >>> <armen.kocharyan at gmail.com> wrote: >>> >>>> Dear Bob, >>>> >>>> Sorry for the asterisks, it wasn't intentional. I just cut and >>>> paste >>>> from my >>>> notebook. >>>> >>>> I understand. But why would Mathematica give different results for >>>> the >>>> following? >>>> >>>> Integrate[h[t], {t, 0, x}] and Integrate[h[x], x] >>>> >>>> >>>> Regards, >>>> Armen >>>> >>>> 2008/5/27 DrMajorBob <drmajorbob at att.net>: >>>> >>>>> The fundamental theorem of calculus frequently does not apply >>>>> when, as >>>>> in >>>>> your example, Integrate is free to return a Complex-valued >>>>> function. >>>>> >>>>> (A better question is why you put in all those extraneous stars to >>>>> inhibit >>>>> copy/pasting the code.) >>>>> >>>>> Despite that, here's a slightly different approach where the FTC >>>>> does >>>>> apply: >>>>> >>>>> q = -1/2; a = 0; b = 3; >>>>> h[x_] = (1 + x^3)^q; >>>>> f[x_] = Integrate[h[t], {t, 0, x}] >>>>> >>>>> If[x >= -1, (x^3)^(1/3) Hypergeometric2F1[1/3, 1/2, 4/3, -x^3], >>>>> Integrate[1/Sqrt[1 + t^3], {t, 0, x}, Assumptions -> x < -1]] >>>>> >>>>> Res1 = N[f[b] - f[a]] >>>>> Res2 = NIntegrate[h[x], {x, a, b}] >>>>> >>>>> 1.65267 >>>>> >>>>> 1.65267 >>>>> >>>>> Contrast that definition of f with yours: >>>>> >>>>> Integrate[h[x], x] >>>>> >>>>> (1/(3^(1/4) Sqrt[1 + x^3]))2 (-1)^( >>>>> 1/6) Sqrt[-(-1)^(1/6) ((-1)^(2/3) + x)] Sqrt[ >>>>> 1 + (-1)^(1/3) x + (-1)^(2/3) x^2] >>>>> EllipticF[ArcSin[Sqrt[-(-1)^(5/6) (1 + x)]/3^(1/4)], (-1)^(1/3)] >>>>> >>>>> If you plot both versions from a to b (use PlotRange->All to be >>>>> sure), >>>>> you'll find that your version is complex on part of the range >>>>> (between >>>>> 2 and >>>>> 3), but mine is real. >>>>> >>>>> Bobby >>>>> >>>>> >>>>> On Mon, 26 May 2008 05:23:35 -0500, Armen Kocharyan < >>>>> armen.kocharyan at gmail.com> wrote: >>>>> >>>>> *q=-1/2; a=0; b=3;* >>>>>> *h[x_]=(1+x^3)^q;* >>>>>> *f[x_]=Integrate[h[x],x];* >>>>>> *Res1=N[f[b]-f[a]]* >>>>>> *Res2=NIntegrate[h[x],{x,a,b}]* >>>>>> >>>>> >>>>> >>>>> >>>>> -- >>>>> DrMajorBob at longhorns.com >>>>> >>> >>> >>> >>> -- >>> >>> DrMajorBob at longhorns.com >>> >> >> >> > > > > -- > DrMajorBob at longhorns.com

**References**:**Re: Integrate vs NIntegrate***From:*DrMajorBob <drmajorbob@att.net>