Re: [問題] 數值積分

看板Mathematica作者 (Amiba Gelos)時間8年前 (2015/04/29 15:49), 編輯推噓4(4014)
留言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

04/29 15:50, , 1F
囧突然發現原本的問題是global maxima
04/29 15:50, 1F

04/29 15:50, , 2F
還以為是global minima哩
04/29 15:50, 2F

04/29 16:44, , 3F
但還是很酷 這題看起來比較像數學問題
04/29 16:44, 3F

04/30 05:53, , 4F
指數函數和Bessel函數裡面的2*Pi/0.7*0.5(^2)怎麼
04/30 05:53, 4F

04/30 05:53, , 5F
改寫掉的?
04/30 05:53, 5F

04/30 05:55, , 6F
問題是最大值一半的2個z值的距離
04/30 05:55, 6F

04/30 05:56, , 7F
所以若能化成近似函數用FindRoot也可以
04/30 05:56, 7F

04/30 05:57, , 8F
但是那個解析解在z=0好像是無窮大
04/30 05:57, 8F

04/30 12:03, , 9F
Bessel裡的常數可以用rescale r來吸收
04/30 12:03, 9F

04/30 12:04, , 10F
rmax從3552/15變成2368Pi/7
04/30 12:04, 10F

04/30 12:05, , 11F
這個近似函數還是只適用於z>r0,所以只知道極大點z<r0
04/30 12:05, 11F

04/30 12:08, , 12F
不過一樣可以對z->0做展開, 但適用範圍會是z<某一定值
04/30 12:08, 12F

04/30 12:10, , 13F
剛剛試了一下對6階展開這個定值大概是z~50@@
04/30 12:10, 13F

04/30 12:11, , 14F
看起來zmax=0 or 50~1000之間
04/30 12:11, 14F

04/30 15:14, , 15F
指數函數改寫的z'=2*Pi/0.7*z*0.5^2 , 這邊的z除以
04/30 15:14, 15F

04/30 15:15, , 16F
2*Pi/0.7*0.5^2 , 就是我的真正z?
04/30 15:15, 16F

04/30 15:23, , 17F
複變的指數函數是週期函數,應該只要用z=0就好了?
04/30 15:23, 17F

04/30 15:24, , 18F
用z=0後的簡化結果 http://goo.gl/4BNrI4
04/30 15:24, 18F
文章代碼(AID): #1LG8piRZ (Mathematica)
文章代碼(AID): #1LG8piRZ (Mathematica)