Re: AiryAi

*To*: mathgroup at smc.vnet.net*Subject*: [mg75692] Re: AiryAi*From*: Roman <rschmied at gmail.com>*Date*: Wed, 9 May 2007 04:15:25 -0400 (EDT)*References*: <f1hlg0$p7q$1@smc.vnet.net><f1mstk$rsi$1@smc.vnet.net>

Great idea to use SequenceLimit! That should be stable, despite of what I said. You can get good estimates for the zero crossings from the asymptotic expansion I gave, z[n_] = -((-1 + 4*n)^(2/3)*(3*Pi)^(2/3))/4 and refine them by numerical search like you do. Cheers! Roman. On May 7, 11:52 am, dimitris <dimmec... at yahoo.com> wrote: > Dear Roman. > > Thanks for all of your information. > I really appreciate your response. > > Finally, trying more, I succeeded in getting a desirable numerical > result. > > In[2]:= > g[o_] := AiryAi[o] > > In[3]:= > dat = Range[-30, -2, 0.5]; > > In[4]:= > zer = Reverse[Append[Union[(o /. FindRoot[g[o] == 0, {o, #1}] & ) /@ > dat, > SameTest -> (#1 - #2 < 10^(-12) & )], 0]] > > Out[4]= > {0,-2.33811,-4.08795,-5.52056,-6.78671,-7.94413,-9.02265,-10.0402,-11.0085, - > \ > 11.936,-12.8288,-13.6915,-14.5278,-15.3408,-16.1327,-16.9056,-17.6613,-18.\ > 4011,-19.1264,-19.8381,-20.5373,-21.2248,-21.9014,-22.5676,-23.2242,-23.871 = > 6,- > \ > 24.5103,-25.1408,-25.7635,-26.3788,-26.987,-27.5884,-28.1833,-28.772,-29.35 = > 48,\ > -29.9318} > > In[4]:= > nint[i_] := -NIntegrate[g[u], {u, zer[[i]], zer[[i + 1]]}, > PrecisionGoal -> 20, WorkingPrecision -> 40] > > In[5]:= > InputForm[SequenceLimit[FoldList[Plus, 0, Table[nint[i], > {i, 1, Length[zer] - 1}]]]] > > Out[5]//InputForm= > 0=2E6666666666666666665047043643`18.362829735502853 > > In[8]:= > Integrate[g[o], {o, -Infinity, 0}] > (N[#1, 20] & )[%] > > Out[8]= > 2/3 > > Out[9]= > 0=2E66666666666666666666666666666666666667`20. > > Dimitris > > =CF/=C7 Roman =DD=E3=F1=E1=F8=E5: > > > Dimitris, > > > The Airy function has terrible convergence: for x -> -Infinity, you > > have, to very good accuracy, > > > AiryAi[x] ~~ Sin[Pi/4 - (2*Sqrt[-x]*x)/3]/(Sqrt[Pi]*(-x)^(1/4)) > > > Just plot the two; they differ by very little for x < -20 or so. Let's > > call f[x] the approximation. > > > So you have an oscillatory function decreasing as slowly as (- > > x)^(-1/4). Numerically integrating this requires adding an infinity of > > terms of alternating signs (the integrals over the positive and > > negative sine arcs, respectively), which have to very delicately > > cancel out in order to produce a small result. The reason why this > > cancellation only works analytically but not numerically in this case > > is that while f[x] is integrable over (-Infinity,0], its square f[x]^2 > > is not. This is like the classical calculus example of a series that > > is convergent but not absolutely convergent, giving rise to numerical > > hell. The sum of all the positive terms is +Infinity, and the sum of > > all negative terms is -Infinity, but the total sum of all terms (the > > integral) is 2/3 (or Sqrt[2/3] for f[x]). Maybe you remember from > > calculus class that by re-ordering the terms in such a series you can > > get any final result you desire, i.e., the result of numerically > > integrating will depend very much on how you integrate (which is bad), > > so I'm not sure that in practice you could ever get a numerical > > integral of your Airy function without analytical insights. > > > Another way of seeing the problem is that as the numerical integration > > progresses from 0 backward to -Infinity, the resulting integral may > > well converge in theory, but the estimated numerical error increases > > the longer we integrate, instead of decreasing or stabilizing as is > > required for a well-determined result. > > > Roman. > > > On May 5, 12:15 pm, dimitris <dimmec... at yahoo.com> wrote: > > > Hello. > > > > Consider the integral > > > > In[678]:= > > > f = HoldForm[Integrate[AiryAi[o], {o, a, b}]] > > > > Then > > > > In[679]:= > > > {ReleaseHold[f /. {a -> 0, b -> Infinity}], (Rationalize[#1, 10^(-12)] > > > & )[ > > > ReleaseHold[f /. Integrate -> NIntegrate /. {a -> 0, b -> > > > Infinity}]]} > > > > Out[679]= > > > {1/3, 1/3} > > > > Hence, the symbolic and the numerical result agree. > > > > Next, > > > > In[682]:= > > > ReleaseHold[f /. {a -> -Infinity, b -> 0}] > > > N@% > > > > Out[682]= > > > 2/3 > > > Out[683]= > > > 0.6666666666666666 > > > > I think I can trust the last analytic result. > > > > However, no matter what options I set for NIntegrate, I could get a > > > satisfactory > > > result from its application for the integral in (-infinity,0]. > > > > I really appreciate any help. > > > > Thanks > > > Dimitris