仔细想想其实这个问题是有解析解的:)
1. z控制了被积函数的震荡频率, 如果震荡频率远快于Bessel的震荡频率的话, 由RPA可知
积分值会趋近于0
2. 因为此函数为正定函数,故可知极值会位在无穷远上
我们可以用MMA来验证一下这个论述是否正确
首先利用BesselJ[0,z] = Integrate[Exp[-I z Sin[t]],{t,-Pi,Pi}]/(2Pi)将原积分式
Integrate[r Abs[Integrate[Exp[I z p^2]BesselJ[0,r p] p,{p,0,1}]]^2
,{r,0,2368 Pi/7}]
换成
Integrate[(R^2+I^2)r,{r,0,2368 pi/7}]
R = Integrate[Cos[z p^2]Exp[-I r p Sin[t]] p,{t,-Pi,Pi},{p,0,1}]/(2Pi)
I = Integrate[Sin[z p^2]Exp[-I r p Sin[t]] p,{t,-Pi,Pi},{p,0,1}]/(2Pi)
这样就可以先对p做积分得
R = Integrate[IntR,{t,-Pi,Pi}]
I = Integrate[IntI,{t,-Pi,Pi}]
IntR = r y ( (FresnelS[r y (2Pi z)^(-1/2)] - FresnelS[(r y+2z)(2Pi z)^(-1/2)])
Sin[r^2 y^2/(4z)] + Cos[r^2 y^2/(4z)] (FresnelC[r y (2Pi z)^(-1/2)]
- FresnelC[(r y+2z)(2Pi z)^(-1/2)]) ) (32Pi z^3)^(-1/2)
+ Exp[-I r y] Sin[z] / (4Pi z)
IntI = r y ( (FresnelS[r y (2Pi z)^(-1/2)] - FresnelS[(r y+2z)(2Pi z)^(-1/2)])
Cos[r^2 y^2/(4z)] - Sin[r^2 y^2/(4z)] (FresnelC[r y (2Pi z)^(-1/2)]
- FresnelC[(r y+2z)(2Pi z)^(-1/2)[) ) (32Pi z^3)^(-1/2)
+ (1-Exp[-I r y] Cos[z])/ (4Pi z)
y = Sin[t]
我这里有做点手动的简化, 如果直接用MMA跑的话会出来error functions (with complex
argument), 看起来很丑Orz (实际上这个简化不会影响对t积分完后的结果)
然后取z->inf的泰勒展开到第2阶得
IntR ~ Exp[-I r y] Sin[z] / (4Pi z) + r y ( 2r y - (2Pi z)^(1/2)
- 2 Sin[r y +z] ) / (16Pi z^2)
IntI ~ (1-Exp[-I r y] Cos[z])/ (4Pi z) + r y ( - (2Pi z)^(1/2)
+ 2 Cos[r y +z] ) / (16Pi z^2)
现在被积函数终于简化到可以对t积分了
R ~ BesselJ[0,r] Sin[z] /(2z) + ( r^2 - 2r BesselJ[1,r] Cos[z] ) /(8z^2)
I ~ (1-BesselJ[0,r] Cos[z])/(2z) + ( - 2r BesselJ[1,r] Sin[z] ) /(8z^2)
最后取平方和对r积分到r0得
r0( r0(BesselJ[0,r0]^2 + BesselJ[1,r0]^2) - 4BesselJ[1,r0] Cos[r0]) /(8z^2)
+ O(r0^3 / z^3)
公式适用范围为 z>>1 and z>>r0
显然当z趋近于无限大, 整个函数趋近于0, 故极值位在无穷远处
与数值结果相比较
http://i.imgur.com/5Phg5cJ.png
可以发现 1.其实z>=r0就已经很准了 2.适用范围其实是 (z>>1 and z>=r0) or z>>r0
如果补上O(r0^3 / z^3)项的话适用范围可以一路拓展到近似公式的first maximam
最后附上MMA code
公式推导
zSeries[F_, n_] := Normal[ Series[ F ,{z, \[Infinity], n} ] ]
dif[F_] := (F /. p -> 1) - (F /. p -> 0)
Int[G_] := ( Integrate[ zSeries[ dif[ # - (# /. Erfi[F_] -> 0) ], 2 ]
// FullSimplify , {y, -1, 1} ]
+ Integrate[ dif[ # /. Erfi[F_] -> 0 ], {y, -1, 1} ] ) & [
Integrate[ G D[ ArcSin[y], y ] / \[Pi], p ] // ExpandAll ]
lim = zSeries[ zSeries[ Integrate[ zSeries[
Int[ Cos[z p^2] E^(-I r p y) p ]^2 + Int[ Sin[z p^2] E^(-I r p y) p ]^2
, 2 ] r, {r, 0, r0} ], 4 ], 2 ]
画图(Warning: q的最大值会严重影响速度, 建议先从5开始试试看)
zmax = 1000; pt = 7;
r0list = Table[ RandomReal[ {0.99, 1.01} ] 2^q, {q, 2, 5} ];
dolist = Flatten[ Table[ {r0, r0^i}, {r0, r0list}, {i, {3/4} ~Join~
Range[ 1, Log[r0, zmax], (Log[r0, zmax] - 1)/(pt - 2) ]} ], 1 ] // N;
NInt[F_, x_, x0_, n_] := NIntegrate[ F, {x, 0, x0}, PrecisionGoal -> n,
AccuracyGoal -> n, Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}]
AbsoluteTiming[ ParallelTable[ { dl[[2]], Quiet[ NInt[ Abs[ NInt[
E^(I dl[[2]] p^2) BesselJ[0, r p] p, p, 1, 7 ] ]^2 r, r, dl[[1]], 5 ] ] }
, {dl, dolist} ] ];
Print["Took ", %[[1]], " sec on a ", $KernelCount, "-Thread Machine"];
Show[ LogLogPlot[ #, {z, dolist[[1, 2]], zmax}, PlotRange -> {10^-6, 1},
AxesLabel -> {z, None}, PlotLegends -> ( "\!\(\*SubscriptBox[ \(r\), \(0\) ]
\)\[Rule]" ~~ ToString[#] & /@ r0list ) ] & [ lim /. r0 -> r0list ],
ParametricPlot[ {Log[zmax], f}, {f, Log[10^-6], 0}, PlotStyle -> LightGray ],
LogLogPlot[ lim /. r0 -> z, {z, dolist[[1, 2]], 1.2 zmax},
PlotStyle -> Black, PlotRange -> All,
PlotLegends -> {"\!\(\*SubscriptBox[ \(r\), \(0\) ]\)\[Rule]z"} ],
ListLogLogPlot[ Partition[ %%[[2]], pt ], PlotMarkers -> Automatic ] ]
要给出一个rigorous的证明其实不容易, 不过原问题r0>>1, 用3阶近似解的话第一极大点
的z远大于1, 所以没意外的话大概就不会有root了