MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: AiryAi

  • To: mathgroup at smc.vnet.net
  • Subject: [mg75826] Re: AiryAi
  • From: m.r at inbox.ru
  • Date: Fri, 11 May 2007 05:32:52 -0400 (EDT)
  • References: <f1hlg0$p7q$1@smc.vnet.net><f1mstk$rsi$1@smc.vnet.net>

In version 6 you can use AiryAiZero:

In[1]:= ii = NIntegrate[AiryAi[z], {z, AiryAiZero[1], 0}] +
  (NIntegrate[AiryAi[z], {z, #[[2]], #[[1]]}]& /@
    Partition[AiryAiZero[Range@ 100], 2, 1] //
      Accumulate);

In[2]:= SequenceLimit[ii]

Out[2]= 0.66666667

Another way is to convert AiryAi to BesselJ:

In[3]:= ee = 1/3 Sqrt[-z] (BesselJ[-1/3, 2/3 (-z)^(3/2)] +
  BesselJ[1/3, 2/3 (-z)^(3/2)]);

In[4]:= FullSimplify[AiryAi[z] - ee, z < 0]

Out[4]= 0

In[5]:= ee2 = (ee /. z-> -t^(2/3)) D[-t^(2/3), t]

Out[5]= -2/9 (BesselJ[-1/3, (2 t)/3] + BesselJ[1/3, (2 t)/3])

In[6]:= NIntegrate[-ee2, {t, 0, Infinity},
  Method -> {"OscillatorySelection", "TermwiseOscillatory" -> True}]

Out[6]= 0.66666667

Maxim Rytin
m=2Er at inbox.ru

On May 7, 4:52 am, dimitris <dimmechan 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.008=
5,=AD-
> \
> 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.8=
71=AD=
> 6,-
> \
> 24.5103,-25.1408,-25.7635,-26.3788,-26.987,-27.5884,-28.1833,-28.772,-29.=
35=AD=
> 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,
>
> > TheAiryfunction 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 yourAiryfunction 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




  • Prev by Date: Re: Mathematica v6: Slower in the following fields
  • Next by Date: Solve & RotationMatrix
  • Previous by thread: Re: AiryAi
  • Next by thread: Re: AiryAi