FC2ブログ

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
スポンサー広告 | --/--/-- --:--

3次ルンゲ・クッタ法

お次は、3次ルンゲ・クッタ

とても挫折しそうな気分なのですが、
前向きな姿勢を込めて、とりあへず、タイトルだけ入れて置きます(笑)

まあ、たまには筋トレもいいかも・・・orz

今度は、以下のような形を仮定する。\[
y_{n+1} = y_n + h(w_1 k_1 + w_2 k_2 + w_3 k_3)
\tag{1}
\]ここで、\[
k_1 = f(x_n, y_n)
\tag{2.1}
\]\[
k_2 = f(x_n + a_2 h, y_n + b_{21} k_1 h)
\tag{2.2}
\]\[
k_3 = f(x_n + a_3 h, y_n + b_{31} k_1 h + b_{32} k_2 h)
\tag{2.3}
\]
$k_2$ をテーラー展開。
(1) の式で、h の3次まで欲しいわけだから、h の2次まで。\[
k_2 = f + (a_2 f_x + b_{21} k_1 f_y)h \\
+ \left( \frac{a_2^2}{2} f_{xx} + a_2 b_{21} k_1 f_{xy} + \frac{b_{21}^2 k_1^2}{2} f_{yy} \right) h^2
+ \mathcal{O}(h^3)
\tag{3}
\]
2変数のテーラー展開なんて、普段やらないから、合ってるかどうか怪しいですが・・・^^;

さて、$k_3$ についても同様に進める(めんどくさそうだが・・・)\[
k_3 = f + \{ a_3 f_x +( b_{31} k_1 + b_{32} k_2 ) f_y \} h \\
+ \left[ \frac{a_3^2}{2} f_{xx} + a_3 ( b_{31} k_1 + b_{32} k_2 ) f_{xy} + \frac{1}{2}
( b_{31} k_1 + b_{32} k_2 )^2 f_{yy} \right] h^2 \\
+ \mathcal{O}(h^3)
\tag{4}
\]
(3) を (4) に代入して、h の2次の項までを考えると、\[
k_3 = f + \{ a_3 f_x +( b_{31} + b_{32} ) f f_y \} h \\
+ \left[ \frac{a_3^2}{2} f_{xx} + a_3 ( b_{31} + b_{32} ) f f_{xy} + \frac{1}{2}
( b_{31} + b_{32} )^2 f^2 f_{yy} \right. \\
\left. + b_{32} (a_2 f_x + b_{21} f f_y) f_y
\right] h^2 \\
+ \mathcal{O}(h^3)
\tag{5}
\]
(3) と (5) を (1) に入れて整理する。
まず、h の1次の項の係数は、\[
(w_1 + w_2 + w_3) f
\tag{6}
\]2次の項の係数は、\[
w_2 (a_2 f_x + b_{21} f f_y)
+ w_3 \{ a_3 f_x +( b_{31} + b_{32} ) f f_y \}
\tag{7}
\]3次の項の係数は、\[
w_2 \left( \frac{a_2^2}{2} f_{xx} + a_2 b_{21} f f_{xy} + \frac{b_{21}^2 f^2}{2} f_{yy} \right) \\
+ w_3 \left[ \frac{a_3^2}{2} f_{xx} + a_3 ( b_{31} + b_{32} ) f f_{xy} + \frac{1}{2}
( b_{31} + b_{32} )^2 f^2 f_{yy} \right. \\
\left. + b_{32} (a_2 f_x + b_{21} f f_y) f_y
\right]
\tag{8}
\]

一方、真の解 $y(x)$ のテーラー展開を3次の項まで書くと、\[
y_{n+1} = y_n + y' h + \frac{1}{2} y'' h^2 + \frac{1}{6} y''' h^3 + \mathcal{O}(h^4)
\tag{9}
\]微分方程式から、\[
y' = f
\tag{10.1}
\]\[
y'' = f_x + f_y y' = f_x + f_y f
\tag{10.2}
\]\[
y''' = f_{xx} + f_{xy} y' + f_{yx} f + f_{yy} y' f + f_y f_x + f_y f_y y' \\
= f_{xx} + 2f_{xy} f + f_{yy} f^2 + f_x f_y + f_y^2 f
\tag{10.3}
\]
(1) と (9) が h の3次の項まで一致するためには、
(6) ~ (8) の係数と (10.1)~(10.3) の係数が一致すればよい。
・・・というわけで、一個一個比較していくと、
1次の項から、\[
w_1 + w_2 + w_3 = 1
\tag{11.1}
\]2次の項から、\[
w_2 a_2 + w_3 a_3 = \frac{1}{2}
\tag{11.2}
\]\[
w_2 b_{21} + w_3(b_{31} + b_{32}) = \frac{1}{2}
\tag{11.3}
\]3次の項から、\[
w_2 a_2^2 + w_3 a_3^2 = \frac{1}{3}
\tag{11.4}
\]\[
w_2 a_2 b_{21} + w_3 a_3 (b_{31} + b_{32}) = \frac{1}{3}
\tag{11.5}
\]\[
w_2 b_{21}^2 + w_3 (b_{31} + b_{32})^2 = \frac{1}{3}
\tag{11.6}
\]\[
w_3 a_2 b_{32} = \frac{1}{6}
\tag{11.7}
\]\[
w_3 b_{21} b_{32} = \frac{1}{6}
\tag{11.8}
\]
ふー、疲れた!どうやら、この中から意味があるのは6個らしい・・・

えーと、まず、(11.7) と (11.8) から\[
a_2 = b_{21}
\tag{12}
\]これを使って、(11.5) と (11.6) から\[
a_3 = b_{31} + b_{32}
\tag{13}
\]すると、(11.4) は自動的に成立するから不要になる。
また、(12) と (13) から、(11.2) と (11.3) の一方は自動的に成立。
以上より、結局、以下の6つの式に集約される。
\[
w_1 + w_2 + w_3 = 1
\tag{14.1}
\]\[
a_2 = b_{21}
\tag{14.2}
\]\[
a_3 = b_{31} + b_{32}
\tag{14.3}
\]\[
w_2 a_2 + w_3 a_3 = \frac{1}{2}
\tag{14.4}
\]\[
w_2 a_2^2 + w_3 a_3^2 = \frac{1}{3}
\tag{14.5}
\]\[
w_3 a_2 b_{32} = \frac{1}{6}
\tag{14.6}
\]
未知数8個に対して、条件は6個なので、やはり解は無数に存在する。

ホインの3次公式\[
w_1 = 1/4, \ w_2 = 0, \ w_3 = 3/4, \\
a_2 = 1/3, \ a_3 = 2/3, \\
b_{21} = 1/3, \ b_{31} = 0, \ b_{32} = 2/3
\] として、\[
y_{n+1} = y_n + \frac{h}{4} (k_1 + 3k_3) \\
k_1 = f(x_n, y_n) \\
k_2 = f \left( x_n + \frac{1}{3}h, y_n + \frac{1}{3} k_1 h \right) \\
k_3 = f \left( x_n + \frac{2}{3}h, y_n + \frac{2}{3} k_2 h \right)
\tag{15}
\]

クッタの3次公式\[
w_1 = 1/6, \ w_2 = 2/3, \ w_3 = 1/6, \\
a_2 = 1/2, \ a_3 = 1, \\
b_{21} = 1/2, \ b_{31} = -1, \ b_{32} = 2
\]として、\[
y_{n+1} = y_n + \frac{h}{6} (k_1 + 4k_2 + k_3) \\
k_1 = f(x_n, y_n) \\
k_2 = f \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_1 h \right) \\
k_3 = f \left( x_n + h, y_n - h k_1 + 2 k_2 h \right)
\tag{16}
\]

この調子で行けば、4次は相当大変そうですね・・・

参考文献
[1] 川上一郎 「数値計算の基礎」 平成21年, http://www7.ocn.ne.jp/~kawa1/
[2] D. Young and R. T. Gregory, "A Survey of Numerical Mathematics Vol. I"
スポンサーサイト
ジャンル:[学問・文化・芸術]  テーマ:[数学
数値計算>微分方程式 | コメント(2) | 2014/10/28 20:17
コメント
管理人のみ閲覧できます
このコメントは管理人のみ閲覧できます
鍵コメさんへ
コメントありがとうございます。

(10.2)については、
y' = f(x,y)= f(x, y(x)) ですから、
y"は f を x で偏微分したもの fx と
f を y で偏微分して(fy)、さらに y を x で微分したもの y'
の和になります。

よって、
y" = fx + fy y' となります。

管理者のみに表示

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。