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