FC2ブログ

スポンサーサイト

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

4次ルンゲ・クッタ法 (1) 失敗編

(この記事は書きかけです)
注:書きながら計算しているので、後で大幅に修正する可能性があります!

さて、いよいよ最後の砦!!!
4次ルンゲ・クッタに登頂を試みるとしますか!(大げさ・・・笑)

完遂できる自信はまったくありませんので、途中で挫折する可能性も大です・・・

まずは、以下の形を仮定。\[
y_{n+1} = y_n + h(w_1 k_1 + w_2 k_2 + w_3 k_3 + w_4 k_4)
\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_4 = f(x_n + a_4 h, y_n + b_{41} k_1 h + b_{42} k_2 h + b_{43} k_3 h)
\tag{2.4}
\]
ああ・・・、これは気が遠くなりますな。

文献 [2] にも、3次まで求めたうえで、
Runge-Kutta formulas of higher order can be derived by the above methods.
However, the labor involved and the complexity increases very rapidly.
と書かれているぐらいです^^;
まあ、ここは、ルンゲさんとクッタさんを信用して、
自分では計算しないのが賢明だとは思いますが・・・(笑)

まずは、h の4次まで欲しいんだから、k2 ~ k4 を3次までテーラー展開する。
まずは、k2。\[
k_2 = f + h(a_2 f_x + b_{21} f f_y ) \\
+ \frac{h^2}{2} \left( a_2^2 f_{xx} + 2 a_2 b_{21} f f_{xy} + b_{21}^2 f^2 f_{yy} \right) \\
+ \frac{h^3}{6} \left[ a_2^3 f_{xxx} + 3 a_2^2 b_{21} f f_{xxy} + 3 a_2 b_{21}^2 f^2 f_{xyy} + b_{21}^3 f^3 f_{yyy} \right] \\
+ \mathcal{O}(h^4)
\tag{3}
\] 続いて、k3 を。\[
k_3 = f + h \left[ a_3 f_x + (b_{31} f + b_{32} k_2) f_y \right] \\
+ \frac{h^2}{2} \left[ a_3^2 f_{xx} + 2 a_3 (b_{31} f + b_{32} k_2) f_{xy} + (b_{31} f + b_{32} k_2)^2 f_{yy} \right] \\
+ \frac{h^3}{6} \left[ a_3^3 f_{xxx} + 3 a_3^2 (b_{31} f + b_{32} k_2) f_{xxy} \\
+ 3 a_3 (b_{31} f + b_{32} k_2)^2 f_{xyy} + (b_{31} f + b_{32} k_2)^3 f_{yyy} \right] \\
+ \mathcal{O}(h^4)
\tag{4}
\]
次に、(3) の k2 を (4) に代入して、それぞれの次数に分ける。
さすがに、全部一度に書けないので、次数ごとに書いていくことにする。
まず、1次の項は、k2 の0次項のみが寄与する。
\[
h \left[ a_3 f_x + (b_{31} + b_{32} ) f f_y \right]
\]
2次の項は、元の2次の項に k2 の0次が入ったものと、
元の1次の項に k2 の1次の項が入ったものが寄与。
\[
\frac{h^2}{2} \left[ a_3^2 f_{xx} + 2 a_3 (b_{31} + b_{32} ) f f_{xy} + (b_{31} + b_{32} )^2 f^2 f_{yy} \\
+ 2b_{32} ( a_2 f_x + b_{21} f f_y ) f_y \right]
\]
で、問題は3次の項。
(a) 元の3次項に k2 の0次項が入ったものと、
(b) 元の2次項に k2 の1次項が入ったものと、
(c) 元の1次項に k2 の2次項が入ったものが寄与するが、

(b) においては、元の2次項には k2 の2乗が含まれるので、
ここに、k2の0次項と1次項の交差項による1次項が入ってることを考慮せねばならない
(ややこしい・・・>_<)

注:上記で「○次項」と言ってるのはすべて、
テーラー展開に現れる $h, h^2, h^3, \cdots$ の意味であって、
$k_2, k_2^2, k_2^3, \cdots$ などの意味ではない。

まず、(a)。
\[
\frac{h^3}{6} \left[
a_3^3 f_{xxx} + 3 a_3^2 (b_{31} + b_{32} ) f f_{xxy} \\
+ 3 a_3 (b_{31} + b_{32} )^2 f^2 f_{xyy} + (b_{31} + b_{32} )^3 f^3 f_{yyy}
\right]
\]
次に、そのややこしい (b)。
まず、k2 の1乗の部分のみ。
\[
h^3 \left[
a_3 b_{32} (a_2 f_x + b_{21} f f_y) f_{xy} \\
+ b_{31} b_{32} (a_2 f_x + b_{21} f f_y) f f_{yy}
\right]
\]
k2 の2乗の部分。
\[
h^3 \left[
b_{32}^2 (a_2 f_x + b_{21} f f_y) f f_{yy}
\right]
\]
最後に、(c)。
\[
\frac{h^3}{2} \left[
b_{32} ( a_2^2 f_{xx} + 2 a_2 b_{21} f f_{xy} + b_{21}^2 f^2 f_{yy} ) f_y
\right]
\]
以上をまとめると・・・
\[
k_3 = f + h \left[ a_3 f_x + (b_{31} + b_{32} ) f f_y \right] \\
+ \frac{h^2}{2} \left[
a_3^2 f_{xx} + 2 a_3 (b_{31} + b_{32} ) f f_{xy} \\
+ (b_{31} + b_{32} )^2 f^2 f_{yy} \\
+ 2b_{32} ( a_2 f_x + b_{21} f f_y ) f_y \right] \\
+ \frac{h^3}{6} \left[
a_3^3 f_{xxx} + 3 a_3^2 (b_{31} + b_{32} ) f f_{xxy} \\
+ 3 a_3 (b_{31} + b_{32} )^2 f^2 f_{xyy} \\
+ (b_{31} + b_{32} )^3 f^3 f_{yyy} \\
+ 6 a_3 b_{32} (a_2 f_x + b_{21} f f_y) f_{xy} \\
+ 6 b_{31} b_{32} (a_2 f_x + b_{21} f f_y) f f_{yy} \\
+ 6 b_{32}^2 (a_2 f_x + b_{21} f f_y) f f_{yy} \\
+ 3 b_{32} ( a_2^2 f_{xx} + 2 a_2 b_{21} f f_{xy} + b_{21}^2 f^2 f_{yy} ) f_y
\right] \\
+ \mathcal{O}(h^4)
\tag{5}
\]
これでようやく、k3 の展開が完了!
この時点で、計算が合ってるとは、とても思えないのですが・・・^^;

次は、k4 のテーラー展開に進みますが、ますます複雑になりそうなので、
ひとまず、ここで記事を区切りたいと思います。

計算が最後まで終わるまでは、合ってるかどうか分からないので、
トップの注意書きはそのまま残しておきます。


参考文献
[1] 川上一郎 「数値計算の基礎」 平成21年, http://www7.ocn.ne.jp/~kawa1/
[2] D. Young and R. T. Gregory, "A Survey of Numerical Mathematics Vol. I"
スポンサーサイト
数値計算>微分方程式 | コメント(0) | 2014/10/30 12:19
コメント

管理者のみに表示

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