作者AmibaGelos (Amiba Gelos)
看板Mathematica
標題Re: [問題] 數值積分
時間Wed Apr 29 15:48:57 2015
仔細想想其實這個問題是有解析解的:)
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了
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 118.166.147.218
※ 文章網址: https://www.ptt.cc/bbs/Mathematica/M.1430293740.A.6E3.html
→ AmibaGelos: 囧突然發現原本的問題是global maxima 04/29 15:50
→ AmibaGelos: 還以為是global minima哩 04/29 15:50
推 willreturn: 但還是很酷 這題看起來比較像數學問題 04/29 16:44
推 obelisk0114: 指數函數和Bessel函數裡面的2*Pi/0.7*0.5(^2)怎麼 04/30 05:53
→ obelisk0114: 改寫掉的? 04/30 05:53
→ obelisk0114: 問題是最大值一半的2個z值的距離 04/30 05:55
→ obelisk0114: 所以若能化成近似函數用FindRoot也可以 04/30 05:56
→ obelisk0114: 但是那個解析解在z=0好像是無窮大 04/30 05:57
→ AmibaGelos: Bessel裡的常數可以用rescale r來吸收 04/30 12:03
→ AmibaGelos: rmax從3552/15變成2368Pi/7 04/30 12:04
→ AmibaGelos: 這個近似函數還是只適用於z>r0,所以只知道極大點z<r0 04/30 12:05
→ AmibaGelos: 不過一樣可以對z->0做展開, 但適用範圍會是z<某一定值 04/30 12:08
→ AmibaGelos: 剛剛試了一下對6階展開這個定值大概是z~50@@ 04/30 12:10
→ AmibaGelos: 看起來zmax=0 or 50~1000之間 04/30 12:11
推 obelisk0114: 指數函數改寫的z'=2*Pi/0.7*z*0.5^2 , 這邊的z除以 04/30 15:14
→ obelisk0114: 2*Pi/0.7*0.5^2 , 就是我的真正z? 04/30 15:15
推 obelisk0114: 複變的指數函數是週期函數,應該只要用z=0就好了? 04/30 15:23