oscillatory integrals
- To: mathgroup at smc.vnet.net
- Subject: [mg69451] oscillatory integrals
- From: dimmechan at yahoo.com
- Date: Wed, 13 Sep 2006 04:01:17 -0400 (EDT)
Hello to all.
***Let me consider the integral of the function q(x,r) (see below) over
the range {x,0,Infinity}, for various values of r.
In[17]:=$Version
Out[17]=5.2 for Microsoft Windows (June 20, 2005)
In[18]:=Clear["Global`*"]
***The function q(x,r) is defined as follows:
In[19]:=q[r_,x_]:=(f[x]/g[x])*BesselJ[0,r*x]
***where
In[20]:=f[x_]:=x*Sqrt[x^2+1/3]*(2*x^2*Exp[(-1/5)*Sqrt[x^2+1]]-(2x^2+1)Exp[(-1/5)*Sqrt[(x^2+1/3)]])
In[21]:=g[x_]:=(2x^2+1)^2-4x^2*Sqrt[x^2+1/3]Sqrt[x^2+1]
***The case r=2 was first considered by Longman on a well celebrated
paper (Longman 1956). I will consider first this case also.
***From the following plot one can see that the integrand is an
oscillatory function convergent to zero for large arguments.
In[22]:=Plot[q[2,x],{x,0,20},PlotPoints\[Rule]1000]
***With its default settings NIntegrate fails to give any result
In[23]:=NIntegrate[q[2,x],{x,0,Infinity}]
NIntegrate::singd: NIntegrate's singularity handling has failed at
point {x}={4.33675*10^14}
for the specified precision goal. Try using larger values for any of
$MaxExtraPrecision or the options WorkingPrecision, or SingularityDepth
and MaxRecursion. More...
NIntegrate::inum : (***I drop the rest of the message***)
Out[23]=NIntegrate[q[2,x],{x,0,8}]
***However increasing the SingularityDepth (the default is 4) gives a
result
In[24]:=NIntegrate[q[2,x],{x,0,Infinity},SingularityDepth\[Rule]6]
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy
after \
7 recursive bisections in x near x = 50.2`. More...
Out[24]=-0.0267271
***Although from the following comands it is clear that in latter case
NIntegrate samples more points, I do not understand why with the
default value of the SingularityDepth NIntegrate fails to give an
answer. Note that for the first message above (NIntegrate::singd:)
there is still no notes in the Online Documentation.
In[26]:=Block[{Message},ListPlot[Reap[NIntegrate[q[2,x],{x,0,Infinity},
SingularityDepth\[Rule]6,EvaluationMonitor\[RuleDelayed]Sow[x]]][[2,1]]]]
In[27]:=Block[{Message},ListPlot[Reap[NIntegrate[q[2,x],{x,0,Infinity},
SingularityDepth\[Rule]6,EvaluationMonitor\[RuleDelayed]Sow[x]]][[2,1]]]]
***Increasing also the maximum number of recursive subdivisions gives a
very reliable result.
In[24]:=NIntegrate[q[2,x],{x,0,Infinity},SingularityDepth\[Rule]20,MaxRecursion\[Rule]20,\
WorkingPrecision\[Rule]22]//Timing
Out[24]={1.406 Second,-0.0266089981279}
***while the following commands dealing with the sampled points
In[29]:=
Length[Reap[NIntegrate[
q[2,x],{x,0,Infinity},
SingularityDepth\[Rule]20,MaxRecursion\[Rule]20,WorkingPrecision\
\[Rule]22,EvaluationMonitor\[RuleDelayed]Sow[x]]][[2,1]]]
Take[Sort[Reap[NIntegrate[q[2,x],{x,0,Infinity},SingularityDepth\[Rule]20,\
MaxRecursion\[Rule]20,WorkingPrecision\[Rule]22,
EvaluationMonitor\[RuleDelayed]Sow[x]]][[2,1]]],-10]
Out[29]=2865
Out[30]=
{1979.89661682140809934,2514.12574728867068961,3744.0426990391850147,3787.\
78762150909090052,5029.2514945773413792,7489.085398078370029,10059.\
5029891546827584,14979.170796156740059,29959.341592313480118,59919.\
68318462696024}
***Next the option Method\[Rule]Oscillatory will be employed
In[32]:=NIntegrate[q[2,x],{x,0,8},Method\[Rule]Oscillatory,WorkingPrecision\[Rule]22]//Timing
8::indet: Indeterminate expression 0\8 encountered. More...
8::indet: Indeterminate expression 0\8 encountered. More...
Out[32]={1.187 Second,-0.026608998128}
***I do not understand why exist here the warning messages
(8::indet:).
Note that despite the presence of the message, the result is very
accurate.
***Now I want to plot the function NIntegrate[q[r,x],{x,0,Infinity}] in
the range {r,0,3}. What is the more reliable method to follow to get
what I want?
I simply executed
Plot[NIntegrate[q[r,x],{x,0,8},Method\[Rule]Oscillatory,
WorkingPrecision\[Rule]30],{r,0,3}]
but although I got a plot, I need considerable time and there were a
lot of warning
messages so I believe this is not the case here.
***Next consider the function h[r,x] which is defined as follows:
In[34]:=h[r_,x_]:=(f[x]/g[x])*BesselJ[1,r*x]
***For this function, we have e.g.
In[34]:=NIntegrate[h[2,x],{x,0,Infinity},SingularityDepth\[Rule]20,MaxRecursion\[Rule]20,\
WorkingPrecision\[Rule]22]//Timing
Out[34]={1.172 Second,-0.147430035385}
***and
In[35]:=NIntegrate[h[2, x], {x, 0, 8}, Method ->Oscillatory,
WorkingPrecision -> 22] // Timing
8::indet: Indeterminate expression 0\8 encountered. More...
8::indet: Indeterminate expression 0\8 encountered. More...
Out[35]={1.078 Second, -0.147430035385}
***which clearly shows the reliability of Method ->Oscillatory.
***I want also here the plot of NIntegrate[h[r,x],{x,0,Infinity}] in
the range {r,0,3}. Because of BesselJ[1,0]=0, I am a little worry how I
will treat
the point r=0.
***Any suggestions?
***Thanks in advance for any assistance.
Dimitrios Anagnostou
NTUA