Taylor展開法による偏微分方程式の数値解法

全文

(1)情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-151 No.14 2015/10/1. Taylor 展開法による偏微分方程式の数値解法平. 山. the spatial direction, it is approximated with high accuracy to partial differential equations in a system of ordinary differential equations, It can be solved by A stable Taylor expansion method, we can get highly accurate results. In the diffusion equation of example, we can get good results whose the absolute error is less than or almost equal to 10−10 .. 弘†1. Taylor 級数の四則演算および関数は C++言語によって容易にできる。四則演算、関数、条件文等で記述された C++言語で定義された関数は容易に Taylor 展開できる。解は任意次数まで計算できるので、Runge-Kutta に代わる任意次数の公式として使うことができる。Taylor 級数を使えば、誤差評価も容易に行え、許容誤差内の適切なステップサイズを容易に求められる。さらに、べき級数を Pad´ e 展開に変換し、それを利用すると任意次数でＡ安定な常微分方程式を解く数値計算法を与える。偏微分方程式を空間的に差分化し、得られる連立常微分方程式を時間方向にべき級数法を適用して解くことを提案する。この方法を使うと安定で精度の高い計算ができる。本文では、空間方向に精度の高いコンパクト差分近似法を使って、偏微分方程式を連立常微分方程式で高精度で近似し、それを A 安定な Taylor 展開法で解き、精度の高い計算が出来た。例題の拡散方程式では、絶対誤差が 10−10 程度以下の計算が出来た。. 1. はじめに有限項で打ち切った Taylor 級数の四則演算、Taylor 級数の関数演算は容易に定義することができる。それを使えば、プログラムの形で与えられた任意の関数をわずかな変更で. Taylor 級数を計算するプログラムに変更3) できる。この変更作業量は、単精度のプログラムを倍精度のプログラムに変更する作業量程度である。Taylor 展開法は、数学的には自動微分7) の一種である。ここで扱う Taylor 級数は、浮動小数点を係数とする有限項で打ち切った級数である。偏微分方程式を、空間的に差分化すると、偏微分方程式は多数の未知数を持つ連立常微分方程式で近似できる。この連立常微分方程式の解を Taylor 級数に展開する。さらにこの. Taylor 級数を Padé 展開すると A 安定な計算法2) になる。. Numerical Solution of Partial Differential Equations by the Taylor Series Method Hiroshi Hirayama†1 The arithmetic operations and functions of Taylor series can be facilitated by the C ++ language. The C++ function written by arithmetic operations, functions, and conditional statements, etc. can be easily Taylor expansion. As we can calculate the solution to any degree, we can be used as the higher order formula to replace the Runge-Kutta. Using the Taylor series, error evaluation also easily be carried out, it can be easily obtained the appropriate step size within tolerance. In addition, it should be to convert the series to Pad ’e series, to give a numerical calculation method to solve the A stable ordinary differential equations in any order when you use it. The partial differential equations spatially differencing, we propose to solve by applying the series method to the simultaneous ordinary differential equations in the time direction which is obtained. This method can be highly stable and accurate calculations when you use it. In this paper, by using the compact differential approximation accurate in. ⓒ 2015 Information Processing Society of Japan. 本論文では、常微分方程式の任意次数の A 安定な計算法を空間的に差分化された方程式に適用し、その性質について述べる。各点の解は A 安定な計算なので、発散することはないが、空間方向の分割数を増加させ、精度を上げようとすると、隣接の解が滑らかに接続しなくなる。このような場合、Taylor 展開式の高次の計算が必要であることがわかる。また、元の偏微分方程式を解を精度良く計算するためには、空間方向にも高精度の公式を使う必要がある。例題では、絶対誤差が 10−10 程度以下になる解を計算することも出来た。計算時間は、Intel Core i7-477K 3.5GHz 上で測定した。コンパイラーとしては Microsoft. Visual Stutio 2013 C++を使用した。 cl /EHsc /O. (ファイル名）. として計算した。. †1 神奈川工科大学 Kanagawa Institute of Technology. 1.

(2) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-151 No.14 2015/10/1. を n としたとき、通常 n3 に比例するのに対し、n2 や n log n に比例する手間で計算できる. 2. 偏微分方程式. 高速な計算法が知られている。本論文の計算は、n2 に比例する Levinson-Durbin の方法6). 偏微分方程式は、時間の偏微分項について解かれた形になっている方程式について考え. によって行った。. る。時間の 2 階微分係数を含む場合、高階常微分方程式を連立１階常微分方程式に変形する. 2.3 A 安定性. ように、1 階微分を別の変数に割りあて、2 階微分を割りあてられた変数の 1 階微分係数と. E. Hairer and G. Wanner1) の章にあるように、A 安定性を調べるために使われる微分方. (. すれば可能である。. ). ∂u ∂u ∂ 2 u = f t, x, u, , ,··· (1) ∂t ∂x ∂x2 ここで、u および f は、一般にベクトルである。この偏微分方程式は、空間方向に差分化. 程式の解を (4) のように、分子 M 次、分母 L 次式に Pad´ e 展開したとき、M ≤ L ≤ M + 2 のとき、A 安定であることが知られている。分母の次数が分子の次数と同じか１または２次高い式に変形すれば、A 安定な公式となる。このように変形すれば常微分方程式の初期値問. することによって、時間についての連立常微分方程式になる。以下の式は等間隔の標本点を. 題に対する任意次数の A 安定な公式が得られたことになる。本方法では、連立一次方程式. 使った場合の式である。 ( ) dun (t) un+1 (t) − un (t) un+1 (t) − 2un (t) + un−1 (t) = fn t, x, un (t), , , · · · (2) dt ∆x (∆x)2 ここで、un (t) = u(t, xn )、xn は n 番目の標本点である。この式では、偏微分の差分化と. を解く程度の計算できるため、計算方法としても、単純で高速である。また、常微分方程式. して最も単純な式を使っている。必要ならば、もっと高次の差分式を利用することもできるし、陰的な数値微分式を利用することも可能である。この方程式は、常微分方程式であるから、次の方法で、解 un (t) は時間の Taylor 展開式の形2) で得られる。. 2.1 常微分方程式の Taylor 級数解次のように. dy dx. について解かれた形になっている常微分方程式について考える。 dy = f (x, y) dx ここで、y および f は、一般にベクトルである。この微分方程式の解のべき級数展開は、以下に示す Picard の逐次計算法8) を使うことによって計算する。 ∫ x. y0 = a0 ,. yn = a 0 +. f (x, yn−1 )dx. を解く方法としてよく使われる Runge-Kutta 法は A 安定な計算方法でないことが知られている。. 3. 偏微分方程式の計算数値例として、次のような単純な偏微分方程式を考える。 ∂u 1 ∂2u = ∂t 12 ∂x2. (5). 境界条件：u(t, x = 0) = u(t, x = 1) = 0 初期条件：u(t = 0, x) = sin πx 積分区間：[0, 10] π2. 厳密解：u(t, x) = e− 12 t sin πx. (3). 0. ただし、(3) の計算は x の n 次以上の高次の項は省略して計算すると効率的に計算できる。この省略をおこなっても、最終的計算結果は同じになる。. 2.2 Pad´ e 展開 Taylor 級数は、容易に Padé 展開することができる。Padé 展開とは、以下のようにべき級数展開式を、有理関数に変形したものである。 p0 + p1 x + · · · + pM xM (4) a0 + a1 x + · · · = x 1 + q1 x + · · · + qL L (4) 式の両辺に右辺の分母を掛け、M + L 次の係数まで一致するように、有理関数の係数を決定することによって得られる。このとき、解かなければならない一次方程式の係数は、. この方程式を (2) のように空間方向に差分化すると次のようになる。 dun (t) un+1 − 2un + un−1 = (6) dt 12(∆x)2 また、時間方向にも差分化すると un ((m + 1)∆t) − un (m∆t) un+1 − 2un + un−1 = (7) ∆t 12(∆x)2 4) 差分方程式 (7) を安定的に解けるためには、次の CFL(Courant-Friedrichs-Lewy) 条件を満たさなければならない。. ∆t ≤6 (∆t)2. (8). Toeplitz 行列と呼ばれる特別な行列になる。この場合の一次方程式の計算量は、未知数の数. ⓒ 2015 Information Processing Society of Japan. 2.

(3) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-151 No.14 2015/10/1. もし、空間方向に 20 等分すれば、∆x =. 1 20. であるから、安定的に解くには、 3 = 0.015 200 以下でなければならない。空間の分割数を 2 倍に増加させると安定な時間の間隔は 4 分の. ∆t ≤. 1 0.8 0.6. 1 になるため、計算反復関数が 8 倍になり、計算時間が 8 倍になる。このように安定な陽的計算法では、急激に増加することが知られている。. 0.4 0.2 0 -0.2. 偏微分方程式を近似する連立常微分方程式 (6) は、次のような手順で Taylor 級数解が得 0. られる。初期条件から、n 番目の点における un (t) の定数項がわかる。得られた定数項を (6) の右辺に代入すると、. dun (t) dt. 1. 0.2 2. 3. 0.4 4. 5. の定数項が決定する。それを積分すれば、1 次の係数が決定され. る。さらに得られた 1 次の項の係数を (6) の右辺に代入すると、. dun (t) dt. 6. 0.6 7. 8. 0.8 9. 10 1. の１次の項が決定図 1 20 分割 8 次の計算結果 (分子分母 4 次). する。それを積分すれば、2 次の係数が決定される。これを繰り返すことによって、任意の次数まで Taylor 展開を求めることができる。この Taylor 展開式を Pad´ e 展開し、その式に刻み幅を代入して、次のステップにおける関数値（定数項）を計算する。これの操作を繰り. も、このように簡単に Taylor 展開式を計算するプログラムを作成できる。実際ここで計算. 返せば、偏微分方程式を解くことができる。. したプログラムは、計算効率が良いので、ほとんどこのような方法で作成することになる。. これをプログラムにすると以下のようになる。un (t) を t = t0 において Taylor 展開した. １次元の空間を 20 等分し、∆t = 0.0125 として、常微分方程式の方法で解いてみた。8. とき、i 次の係数を u[n][i] と表現したとき、次のようなプログラムとなる。. 次式で計算した結果を図 1 に示す。A 安定な計算方法なので、一部スパイク状のものが現. for( int k=0 ; k<ORDER ; k++ ) // ORDER 次の係数まで計算. れるが安定に計算できる。この場合、計算は安定であるが、計算精度はあまり良くない。ス. {. パイクが現れないようにするには、計算次数を上げる方法と時間の刻み幅を小さくする方法. u[0][k+1] = 0 ;. // 端点の値. for( int j=1 ; j<NX ; j++ ) // NX は空間方向の分割数. がある。この問題の場合、次数を 10 以上にするとスパイクが現れないことがわかる。その計算結果 (10 次の計算) を図 2 に示す。精度上げるために、分割数を 50 にした。スパイク状の結果が現れなくするには、計算次. {. work[j] = cf*(u[j+1][k]-2*u[j][k]+u[j-1][k] ) ; // cf=1/(12*h*h) 数を 28 次にする必要があった。時間の刻み幅を小さくさせると、計算時間が増加するとい }. う問題が生じる。分割数を 50 にした計算結果を図 5 に示す。20 分割と 50 分割と分割数を. for( int j=1 ; j<NX ; j++ ). 増やし、計算次数を上げても、厳密解との差は 3 桁程度の一致するだけで、精度の良い結果. {. は得られなかった。例題では、区間 [0,10] の計算しか行っていないが、安定な計算では区間. u[j][k+1] = work[j]/(k+1) ; // 積分を行う。 } u[NX][k+1] = 0 ; // もう一方の端点の値 } Taylor 展開式の計算であるが、特別に準備された Taylor 展開のプログラムを使わなくて. ⓒ 2015 Information Processing Society of Japan. [0,100] でも問題なく計算できる。計算は安定で、最もらしい計算結果が得られるが、その精度は 3 桁程度である。さらに精度の良い結果を得るには、偏微分方程式を連立常微分方程式で近似する段階で、さらに近似あげる必要がある。今回はこのために、陰的微分法であるコンパクト差分法5) を利用した。. 3.

(4) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-151 No.14 2015/10/1. 立つ。 (2). (3). fi f h2 + i h3 + O(h4 ) (9) 2! 3! ここで、h は標本点の間隔である。h の代わりに −h を代入すると、次の式が得られる。 (3) (2) f f (1) (10) fi−1 = fi − fi h + i h2 − i h3 + O(h4 ) 2! 3! (9) と (10) を加えると、次の式が得られる。 (6) (4) f f (2) (11) fi+1 + fi−1 = 2fi + fi h2 + 2 i h4 + 2 i h6 + O(h8 ) 4! 6! (4) fi は、fi (2) の 2 階微分であるから、次の陽的な 2 階微分公式が使える。 (2) (2) (2) f − 2fi + fi−1 (4) fi = i+1 + O(h2 ) (12) h2 (12) を (11) に代入すると fi+1 − 2fi + fi−1 1 (2) (2) (2) (f + 10fi + fi+1 ) = + O(h4 ) (13) 12 i−1 h2 各標本点で、(13) の式が成り立つから、3 項だけゼロでない項を持つ連立一次方程式にな (1). fi+1 = fi + fi h +. 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0. 1. 2. 3. 0.6 4. 5. 6. 7. 0.8 8. 9. 10 1. 図 2 20 分割 10 次の計算結果 (分子分母 5 次). る。この方程式の係数行列は、3 重対角行列となる。この行列は、LU 分解しても、非ゼロ係数は増加しないので、簡単に逆行列に相当する式が得られる。これを利用すれば、この方程式は容易に解くことができるため、2 階微分係数. 1 0.8. の計算も容易に高次 (h4 ) で計算できる。. 0.6 0.4. この差分式は、両隣の値を利用して計算している点で、最初の方程式 (6) と同じで扱いや. 0.2 0. すい。(13) の公式を利用すれば、比較的変更も少なく高精度の 2 階微分係数を計算するこ 0.2 0.4 0. 1. 2. 3. とができる。. 0.6 4. 5. 6. 7. 0.8 8. 9. 10 1. 3.2 高次の計算初期条件は、x = 0.5 について、対称性であり特別な解を計算しているように思えるので、初期条件を対称性のない初期条件にして、高精度計算を行った。. 図 3 50 分割 28 次の計算結果 (分子分母 14 次). ここで使用した初期条件は、次のように設定した。初期条件：u(t = 0, x) = sin πx + sin 2πx 積分区間：[0, 10] π2. 厳密解：u(t, x) = e− 12 t sin πx + e−. π2 3. t. sin 2πx. 区間を 500 等分して、28 次の公式を利用して、∆t = 0.004 として計算した。その結果を. 3.1 コンパクト差分法ここで行った計算で使用したコンパクト差分法に説明する。等間隔にある i 番目の標本点 (n). における関数値を fi 、n 階微分係数を fi. と書くと、次の関係式 (Taylor 展開式) が成り. 示す。∆t = 0.005 にすると、スパイク状のものが現れ、精度が著しく損なわれる。それでも見た目は安定で最もらしい計算結果である。計算時間は約 17 秒であった。グラフを書くためファイル出力をすると約 32 秒であった。. ⓒ 2015 Information Processing Society of Japan. 4.

(5) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-151 No.14 2015/10/1. 3.3 A 安定でない計算 Padé 展開式を使えば、A 安定な計算になる。高い精度で計算すれば、安定な計算になるのは当然である。. f (a + h) = f0 + f1 h + f2 h2 + f3 h3 + · · · + fn hn + O(hn+1 ). (14). ここで、最高次項が要求精度 ϵ より小さいとする。. |fn | hn ≤ ϵ. 2 1.5. すなわち. 0. 1 0.5. ϵ (16) |fn | このように、計算ステップサイズを決定すれば、高い精度で計算できる。28 次の Taylor 展. 0 0.4. 0.8 2. 4. 6. 8. n. 開式を利用して、上の例題を計算してみた。この時の ϵ = 10−10 とした。約 70000 ステッ. 0.6 0. √. h≤. 0.2. -0.5. (15). プで同様な計算が出来た。計算時間は、約 38 秒であった。Pad´ e 展開を利用した方法とほ. 1. とんど同じであった。このようにステップサイズをとれば、ほぼ 10 桁程度の結果が得られるはずであるが、(16). 図 4 500 分割 28 次の計算結果 (分子分母 14 次). の条件が Taylor 展開式が精度良く計算できる十分条件ではない。さらに、計算途中の計算結果が最終結果に比べ絶対値があまり大きくないという条件が付く。. 4. 結. 論. A 安定な Taylor 展開法を利用して偏微分方程式を解いた。空間方向は単純な差分によって近似した場合、Taylor 展開の次数が上がると精度の向上が見られるが、空間方向の差分近似が低次のためか、あまり精度良く計算できなかった。空間方向の差分近似を陰的方法を使い、より高次の近似式を使うと絶対誤差で 10−10 の精度で計算ができる。. 2 1.5 1. 0. 0.5 0.2. 0 -0.5. 0.4 0.6 0. 1. 2. 3. 4. 5. また、拡散方程式だけでなく、他の発展型偏微分方程式に対する性能を評価したいと考えている。. 参. 考. 文. 献. 0.8 6. 7. 8. 9. 10 1. 図 5 500 分割 28 次の計算の絶対誤差 (1010 倍). ⓒ 2015 Information Processing Society of Japan. 今後は空間方向の差分近似を高精度化して、この方法を適用する予定である。. 1) Hairer E., Wanner G., Solving Ordinary Differential Equations II, SpringerVerlag,(1991) 2) 平山, 小宮, 佐藤, Taylor 級数法による常微分方程式の解法, 日本応用数理学会論文誌 12(2002), 1–8. 5.

(6) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2015-HPC-151 No.14 2015/10/1. 3) 平山, 舘野, 浅野, 川口, Taylor 級数演算ライブラリの使用法, 東北大学情報シナジーセンター大規模科学計算システム広報 SENAC, 40(2007) 29–68 4) 河村哲也, 流体解析の基礎, 朝倉書店,(2014) 5) 桑原邦郎, 河村哲也, 流体計算と差分法, 朝倉書店,(2014) 6) Press W.H., Flannery B.P., Teukolsky S.A., Vetterling W.T., Numerical Recipes in C, Cambridge University Press, Cambridge(1988) 7) Rall,L. B., Automatic Differentiation Technique and Applications, Lecture Notes in Computer Science, Vol. 120, Springer Verlag, Berlin-Heidelberg-New York, 1981 8) 佐野理, キーポイント微分方程式, 岩波書店, 東京, (1993). ⓒ 2015 Information Processing Society of Japan. 6.

(7)