高振動型関数の数値積分法

全文

(1)情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2013-HPC-139 No.5 2013/5/29. In order to perform this calculation, it is necessary to compute higher order differentiation of the above function. It is shown that the above integration is efficiently calculable using numerical Taylor series in this calculation.. 高振動型関数の数値積分法平. 弘†1. 山. 1. はじめに. 本論文では、Fourier 級数の計算等に現れる次のような形を持つ１次元の高振動関 ∫ b. 数のについて考える。I =. f (x)e. iωg(x). dx. (ω > 0) この積分は ω が大きくなる. a. と、計算が著しく難しくなる。この積分を、1/ω について級数展開すると、ω が大きな値の場合、容易に計算できる。積分区間内で g 0 (x) 6= 0 の場合、g(x) の逆関数が存在するので、その逆関数を利用すれば、上の式は t = g(x) と置けば、次のように変形できる。. ∫. g(b). 1 −1 (g (t))eiωt dt ω > 0 dt g(a) このように変形できれば、部分積分法によって、容易に 1/ω についての級数展開できる。この計算には、関数の高階の微分係数を必要するが、その計算には数値的 Taylor 展開を利用して効率的に計算できることを示した。 I=. f (g −1 (t)). 高振動関数の積分として、次のような数値積を考える。. ∫. b. f (x)eiωx dx. I=. (ω > 0). (1). a. このような積分に関数する数値積分法については、多くの研究2)6)7) が行われている。振動数が小さい場合には、通常の数値積分法で容易に積分可能であるが、振動数が増すと関数計算の回数が増加し計算が著しく難しくなる。以下に積分で ω にいろいろな数値を与えて倍精度で計算した。. ∫. I= 0. 1. sin ωt dx 1 + x2. (2). 数値積分法として二重指数型数値積分公式を利用した。その結果を表 1 に示す。この計算表 1 二重指数型数値積分法による相対誤差と標本点数. Numerical Integration Methods for Highly Oscillatory Function Hiroshi Hirayama†1 In this paper we consider the numerical integration methods for 1-dimensional ∫. ω 10 50 100 500 1000 5000. 積分値 0.14658077366028335000 0.01042276124859460600 0.00571615773800527530 0.00288479728232866230 0.00071839930329983944 0.00018455293532718502. 相対誤差 5.6806e-016 1.59779e-014 4.65837e-014 2.89241e-013 3.04253e-013 1.19047e-011. 標本点数. 103 205 411 1637 3301 4711. 1. f (x)eiωg(x) dx. highly oscillatory integrals of the form I =. (ω > 0) If ω be-. 0. comes large, calculation of above integrations will become remarkably difficult. When the series deployment for above integration of 1/ω is carried out and ω is a big value, it can calculate easily. Since the inverse function of g(x) exists in the integration interval in g 0 (x) 6= 0, if the inverse function is used, the above formula can change as follows, if it places with t = g(x). ∫ g(1) 1 I= f (g −1 (t)) (g −1 (t))eiωt dt ω > 0 dt g(0) if it can change, series of 1/ω can be easily calculated with integration by parts.. ⓒ 2013 Information Processing Society of Japan. 結果から、ω が 1000 を超えれば、標本点数は 3000 を越えかなり計算が難しいことがわかる。また精度も標本点数が増加するとともに悪くなっていくことがわかる。本論文では、ω が大きいとき、効率的に計算する方法を提案する。. †1 神奈川工科大学 Kanagawa Institute of Technology. 1.

(2) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2013-HPC-139 No.5 2013/5/29. ば、そのとき誤差が最小となる。この誤差が十分小さいならば、この級数を積分計算に使用. 2. ω による逆数展開式. できることになる。. (1) の積分は、部分積分法. ∫. f (x)eiωx dx = −i. f (x) iωx i e + ω ω. 実際上は、各項の値の絶対値が増加し始めたら、それまでの計算した値がその級数の値と. ∫. する。誤差は計算した最終項の絶対値程度になる。. f 0 (x)eiωx dx. (3). この公式の実数部分と虚数部分1) をとれば、それぞれ. ∫. を使って、. ∫. b iωx. I=. f (x)e. [. dx. f (x) = −i f (x)eiωx ω. [{. −i. [ = eiωx. ]b. i + ω a. f (x) f 0 (x) + ω ω2 n { ∑ k=0. (−1)n+1 − 2n+2 ω. }. −i(−1)k. ∫. ∫. S=. b. 0. iωx. f (x)e a b. ]. −. eiωx a. 1 ω2. ∫. b. f (x) sin wxdx. (6). a. dx. の積分値が得られる。(6) の C と (6) の S の積分は、一つの積分を計算しても、二つ同時に計算しても計算量はほとんど変わらない。. f 0 (x)eiωx dx. a. f (2k+1) (x) f (2k) (x) + (−1)k 2k+1 ω ω 2k+2. 積分公式の多くの公式は、積分区間 [a, b] の長さに依存する公式が多いが、(4) の公式は、. }]b. 積分区間 [a, b] の長さに依存しない公式となっている。積分区間が短いと高精度で計算が可. (4). 能な公式が多いが、今回の公式は区間の長さに依存しないので区間が大きい問題に有効に働くものと考えられる。逆に積分区間が短い問題は、この公式で計算することは、通常の数値. a b. 積分公式と比較し相対的に得策とはならないように思える。. f (2n+2) (x)eiωx dx a. のように展開できる。この展開式は、漸近展開式であることに注意する必要がある。(4) の最後の積分は、最初の何回かの部分積分法では、徐々に小さくなっていくが、途中から徐々. 3. ω の漸近級数による計算例最初に上げた例 (2) について計算する。. ∫. に大きくなる。(4) の積分項が十分小さくなれば、そのときの漸近級数の値は精度の値を与. 1. I=. える。この性質を使って計算することになる。. 0. ¯ ¯ ∫ ¯ ¯ b−a ¯ ¯ (−1)n+1 b (2n+2) iωx ¯ f (x)e dx¯¯ ≤ n+1 ¯f (2n+2) (ζ)¯ E = ¯ 2n+2 ω ω. sin ωt dx 1 + x2. (7). この積分の f (x) に相当する関数は、何回か微分すると徐々に微分係数が大きくなり発散す. 最終項の絶対値は. ζ ∈ [a, b]. (5). a. もし、f (x) が n 次の多項式ならば、n + 1 回微分すると 0 になるので、多項式のに対しては、厳密な公式による計算になるので、打切誤差が入らない精度の高い計算ができる。単純な指数関数 epx ならば、n 回微分が pn epx になるので、誤差項の大きさは O(( ωp )n ) のオーダーになるので、p < ω ならばいくらでも小さくなる。すなわち漸近級数ではなくなる。単純な分数関数. f (x) cos wxdx. ∫a b. a. =. b. C=. る関数である。このため ω が十分に大きくない場合には、精度良く計算できない場合が考えられる。. f (x) =. 1 1+x2. を x = 0 で Taylor 展開すると次のようになる。計算は 30 次まで計算して. ある。事前に必要な次数がわからないためかなり高い次数まで計算してある。. f(x)=1-x^2+x^4-x^6+x^8-x^10+x^12-x^14+x^16 -x^18+x^20-x^22+x^24-x^26+x^28-x^30 同様に x = 1 で Taylor 展開すると. 1 1+x. である場合、n 回微分すると. (−1)n n! (1+x)n+1. となる。この場合、n! は n. が ω を越えると増加し始め、最終的には発散する。この場合、n = ω あたりまで計算すれ. ⓒ 2013 Information Processing Society of Japan. f(x)= 0.5-0.5*(x-1)+0.25*(x-1)^2-0.125*(x-1)^4+0.125*(x-1)^5 -0.0625*(x-1)^6+0.03125*(x-1)^8-0.03125*(x-1)^9+0.015625*(x-1)^10. 2.

(3) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2013-HPC-139 No.5 2013/5/29. 行っていない。. -0.0078125*(x-1)^12+0.0078125*(x-1)^13-0.00390625*(x-1)^14. 積分区間内で g 0 (x) = 0 になる場所がある場合には、その部分を除いて何個かの区間に分. +0.00195312*(x-1)^16-0.00195312*(x-1)^17+0.000976562*(x-1)^18 -0.000488281*(x-1)^20+0.000488281*(x-1)^21-0.000244141*(x-1)^22. 割し、逆関数による変数変換を行えば、前節までの計算法で計算できる。g 0 (x) = 0 を含む区. +0.00012207*(x-1)^24-0.00012207*(x-1)^25+6.10352e-05*(x-1)^26. 間は、その区間を狭くして、通常の数値積分法で計算すれば、容易に求められる。g 0 (x) = 0. -3.05176e-05*(x-1)^28+3.05176e-05*(x-1)^29-1.52588e-05*(x-1)^30. であることを利用して、鞍点法で計算する方法も研究9) されている。. 係数は 5,6 桁程度の精度で表示されているが、倍精度実数で計算してあるので、計算精度は. t = g(x) とおいて、(8) の式を変数変換すると. ∫. 約 15 桁である。この漸近級数を使って計算すると表 2 のようになる。この展開式は、すべての ω について共通である。. g(b). I=. f (g −1 (x)). g(a). 表 2 ω の漸近級数による計算. ω 10 50 100 500 1000 5000. d −1 (g (x))eiωx dx dx. 相対誤差. 項数. 0.14661439599958586700 0.01042276124859443584 0.00571615773800500382 0.00288479728232782793 0.00071839930329962011 0.00018455293532938208. 4.8760866e-05 1.3238450e-17 8.2908485e-18 3.3507995e-18 1.1315575e-18 2.8335267e-19. 12 20 12 8 6 6. 4.1 逆関数の計算逆関数の Taylor 展開は、直接 Taylor 展開する方法や逆関数が満たす微分方程式の解を. Taylor 展開する方法がある。その他に、Lagrange inversion formula3) がある。この方法によると関数 f (x) の逆関数 g(x) は. g(x) = f (x) − f (0) x g(x). ω = 10 の場合、12 項から絶対値が増加したため、ここで計算を打ち切っている。このと. cn =. きの誤差は倍精度計算としては非常に悪い精度になっている。それ以外は、すべて倍精度の最大精度まで計算が出来ている。ω が 1000 を越えると計算は 6 項程度の計算で済み高精度で計算できることがわかる。この計算を通常の数値積分法で計算すれば、何千個の標本点を計算する必要がある。このような計算では、本方法は、非常に効率的である。. (10). (ω > 0). (8). a 0. ここで、積分区間内で g (x) 6= 0 であるとする。この条件を満たすと関数 g(x) の逆関数が存在する。ここでは実際に、逆関数の Taylor 展開を求め、変数変換を行う。この方法については、Lyness6) などでその考え方およびその方法について述べているが、実際の計算は. ⓒ 2013 Information Processing Society of Japan. (12) x=0. cn (x − f (0))n. (13). n=1. 数 f (x) の逆関数の微分方程式が得られる。. 1 dy = 0 dx f (y). b. dx. (11). 関数 y = f (x) の逆関数は x = f (y) と表すことができる。この式を微分すると、次の関. ここでは、さらに一般的な高振動型数値積分に述べる。ここで考察する積分は、次のよう. f (x)e. ∞ ∑. }¯ ) ¯ d n−1 n ¯ (w(x)) ¯ dx. の数値解法になる。. な形式のものである. I=. f −1 (x) =. 1 n!. {(. として、逆関数を求めることができる。この計算は何次まででも計算できるので、任意次数. 4. 逆関数による変数変換. iωg(x). (9). となる。. 積分値. w(x) =. ∫. (ω > 0). (14). 初期値は、適当に y0 の値を決め、x0 = f (y0 ) として決める。得られた微分方程式は、次の方法で Taylor 展開ができる。一般に常微分方程式は次のように与えられているとする。. dy = f (x, y) dx. y(0) = a0. (15). 3.

(4) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2013-HPC-139 No.5 2013/5/29. このとき、Picard の逐次計算法11) を使って、次のように求められる。. y0 = a 0. ∫. かもしれない。. (16). 逆関数が、解析的方法で陽に求められない場合の例題を上げる。. ∫. x. yn = a 0 +. f (x, yn−1 )dx. (17). 0. x. (22). 0. (17) の反復計算は、通常 1 回の反復計算で解の精度が１次上がるので、この計算を n 回繰. ここで t = g(x) = x + ex と置くと、次のように変形される。. ∫. り返せば通常 n 次の Taylor 展開が得られる。通常解の精度は１次上がるが、2 次上がる例. y(0) = a0. 1+e. eg. I=. もある。次のような方程式の場合、1 回の計算で 2 次の精度が上がる。. dy = xy dx. 1. ex eiω(x+e ) dx. I=. −1. 1. (18). (t). d −1 (g (t))eiωt dt dt. (23). (23) の積分の下限は g(0)、上限は g(1) である。下限は 1 であるから、高振動項部分を除いた被積分関数を x = 1 で Taylor 展開をしなけ. n 回目の反復計算で n 次まで正しい結果が得られるので、Taylor 展開の計算次数を n 次に. ればならない。このためには、g(x) の逆関数を逆関数の展開位置は、元の関数の値になる. して計算する。また、k 回目の反復計算には、k − 1 回目の反復計算と同じ計算が含まれて. ので、x = 0 で g(x) を Taylor 展開し、その逆関数を求めれば、その逆関数は x = g(0) で. いるのでそれを利用すると約 2 倍速さで計算出来る可能性がある。この場合前回の計算の. Taylor 展開される。上の例題では、g(x) を x = 0 で Taylor 展開し、20 次まで表示すると. 途中経過を保存しておく必要があるので、記憶容量を多く必要とする。この方法で計算す. 次のようになる。. れば、逆関数の Taylor 展開も容易にできる。今回の例題は、計算精度が良かった Picard. g(x)= 1+2*x+0.5*x^2+0.166667*x^3+0.0416667*x^4+0.00833333*x^5+0.00138889*x^6. の逐次計算法4) を使って計算を行った。この計算プログラム等は、平山等5) にある。また. +0.000198413*x^7+2.48016e-05*x^8+2.75573e-06*x^9+2.75573e-07*x^10. 8)10). Taylor 展開だけでなく自動微分については、久保田等. を参照してください。. +2.50521e-08*x^11+2.08768e-09*x^12+1.6059e-10*x^13+1.14707e-11*x^14 +7.64716e-13*x^15+4.77948e-14*x^16+2.81146e-15*x^17+1.56192e-16*x^18. 5. 逆関数を利用した数値計算例. +8.22064e-18*x^19+4.11032e-19*x^20. ここでは、(8) で示された次の形の積分を扱う。. ∫. この逆関数を Taylor 展開すると次のようになる。展開位置が x = 1 になっている。. b. f (x)eiωg(x) dx. I=. (ω > 0). (19). a. 次のような例題を考える。. ∫. 1. 2. ex eiω(x+1) dx. I=. (20). これは、解析的に容易に変数変換できて、次のようになる。. I= 1. 4. e. −1+. √. √ x. 2 x. xeiωx dx. -0.000211589*(x-1)^5+3.18739e-05*(x-1)^6+1.76808e-06*(x-1)^7-1.8521e-06*(x-1)^8 +3.53444e-07*(x-1)^9+8.17383e-09*(x-1)^10-2.06245e-08*(x-1)^11. 0. ∫. g^{-1}(x)= 0.5*(x-1)-0.0625*(x-1)^2+0.00520833*(x-1)^3+0.000325521*(x-1)^4. +4.72725e-09*(x-1)^12-5.25484e-11*(x-1)^13-2.54973e-10*(x-1)^14 +6.89708e-11*(x-1)^15-3.00464e-12*(x-1)^16-3.31398e-12*(x-1)^17 +1.05477e-12*(x-1)^18-7.81129e-14*(x-1)^19-4.39179e-14*(x-1)^20 この結果を、f (g −1 (x)) および. (21). d −1 g (x) dx. に代入すると、次の結果が得られる。. result = 0.5+0.125*(x-1)-0.015625*(x-1)^2-0.00130208*(x-1)^3+0.00105794*(x-1)^4. このように逆関数を解析的に陽に書くことが出来る場合、逆関数の Taylor 展開を使わずに. -0.000191243*(x-1)^5-1.23766e-05*(x-1)^6+1.48168e-05*(x-1)^7-3.181e-06*(x-1)^8. 容易に計算できる。逆関数を解析的に求めることが出来る場合でも、逆関数の分岐がたくさ. -8.17383e-08*(x-1)^9+2.2687e-07*(x-1)^10-5.67271e-08*(x-1)^11. んあり、それを考慮しなければならない場合には、単純に Taylor 展開を使った場合が簡単. +6.83129e-10*(x-1)^12+3.56963e-09*(x-1)^13-1.03456e-09*(x-1)^14. である場合がある。(20) の例題では、分岐点は 1 個しかないので解析的に行った法が簡単. +4.80742e-11*(x-1)^15+5.63377e-11*(x-1)^16-1.89859e-11*(x-1)^17. ⓒ 2013 Information Processing Society of Japan. 4.

(5) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2013-HPC-139 No.5 2013/5/29. +1.48415e-12*(x-1)^18+8.78358e-13*(x-1)^19-3.48048e-13*(x-1)^20 このように、最終結果も x = 1 で Taylor 展開されている。同様に x = g(1) = 1 + e でも同様に計算できる。計算結果は 0.00113925295659538916 で、推定誤差= 1.0243692e-17 である。使用した. Taylor 展開の項数=8 であった。8 次までの Taylor 展開があれば計算はできるが、事前に何次まで必要であるかわからないため、現在は 30 次まで Taylor 展開を行っている。この結果は、ほかの計算法とほぼ一致する。同様な方法で、(23) の積分を計算すると、0.00787401314319245878、相対誤差=. 1.1129301e-17、必要な項数=12 であった。上で述べたように、逆関数を指定の位置 s で Taylor 展開するためには、s = g(x) を満たす x を求める必要がある。今回の問題ではこの部分が自明であるから容易に計算できたが、. 12(2002), 1–8 5) 平山, 舘野, 浅野, 川口, Taylor 級数演算ライブラリの使用法, 東北大学情報シナジーセンター大規模科学計算システム広報 SENAC, 40(2007) 29–68 6) James N. Lyness, JamesW. Lottes, Asymptotic expansions for oscillatory integrals using inverse functions, BIT Numer Math 49(2009) 397–417 7) Veerle Ledoux, Marnix Van Daele, Interpolatory Quadrature Rules for Oscillatory Integrals, J. Sci. Comput. 53(2012) 586–607 8) 久保田光一, 伊理正夫, アルゴリズムの自動微分と応用、コロナ社、(1998) 9) Olver, S., Fast, numerically stable computation of oscillatory integrals with stationary points. BIT 50(2010) 149–171 10) Rall,L. B. , Automatic Differentiation-Technique and Applications, Lecture Notes in Computer Science, Vol.120, Springer-Verlag, Berlin-Heidelberg-New York(1981) 11) 佐野理, キーポイント微分方程式, 岩波書店, 東京, (1993). 必要があれば Newton 法か何かで徳必要がある。. 6. 終わりに高振動関数の数値積分法を、理論だけでなく実際に計算を行った。かなり良い結果が得られた。高振動関数の積分について述べられている論文では、微分係数や Taylor 展開を利用して計算してあるものもあるが、微分計算や Taylor 展開は計算機で行うのではなく、手計算で行っているらしくそのため、計算次数が非常に低いものであった。この点は改善できたと考えている。また、これらの論文では、特異点の扱いについては、非常に詳しく書かれているので、それらも数値 Taylor 展開法で計算を行っていきたいと考えている。三角関数の高振動型関数を考えたが、Fouries-Bessel 級数展開に現れる Bessel 関数の高振動型関数の計算も考えて行きたいと考えている。. 参. 考. 文. 献. 1) Abramowitz, M. & Stegun, I.A. Handbook of mathematical functions. Washington, DC: National Bureau of Standards. (1964) 2) Davis, P.J. & Rabinowitz, P. Methods of numerical integration, 2nd edn. Orlando, FL: Academic Press. (1980) 3) Bruijn, N.D.: Asymptotic Methods in Analysis,Dover Pub. (1981). 4) 平山, 小宮, 佐藤, Taylor 級数法による常微分方程式の解法, 日本応用数理学会論文誌,. ⓒ 2013 Information Processing Society of Japan. 5.

(6)