Re: [問題] 數值積分

看板Mathematica作者AmibaGelos (Amiba Gelos)時間9年前 (2015/04/29 15:49)推噓4(4推 0噓 14→)

留言18則, 3人參與討論串2/3 (看更多)

仔細想想其實這個問題是有解析解的:) 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

04/29 15:50, , 1^F

04/29 15:50, 1^F

→

AmibaGelos

04/29 15:50, , 2^F

04/29 15:50, 2^F

推

willreturn

04/29 16:44, , 3^F

04/29 16:44, 3^F