スポンサーサイト

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

二階微分方程式の数値求解 (2)

とりあえず、クラスへの実装の話は置いておいて、
無骨にプログラミングして、計算してみました。

刻みを h = 0.1 に設定して、計算した結果。
偶関数の初期条件の場合。
RungeKuttaCos01.png
奇関数の初期条件の場合。
RungeKuttaSin01.png
赤が計算結果で、青が解析解。
比較してみると、きちんと計算できてそうです。

せっかくなので、刻みをどれだけ粗くできるかチェックしてみることに。

h = 0.2
RungeKuttaCos02.png

h = 0.5
RungeKuttaCos03.png

h = 1
RungeKuttaCos04.png

さすがに、ここまで粗くすると、少しずれが見えてきましたが、
結構、粗くても大丈夫なんですね。
振動周期の1/10程度なら大丈夫ということが分かりました。

というわけで、次回、シュレディンガー方程式の計算を試してみたいと思います。
スポンサーサイト
数値計算>微分方程式 | コメント(0) | 2014/12/17 19:38

二階微分方程式の数値求解 (1)

シュレディンガー方程式を解くために、
まずは、二階の微分方程式を数値的に解く練習。

最も簡単なもので、\[
y'' = -y
\tag{1}
\]を解くことにする。

解析解は、奇関数の初期条件 $y(0) = 0$, $y'(0) = 1$ で\[
y = \sin x
\tag{2.1}
\]偶関数の初期条件 $y(0) = 1$, $y'(0) = 0$ で\[
y = \cos x
\tag{2.2}
\]
$y' = z$ とおくと、(1) は以下の連立微分方程式に帰着。\[
y' = z
\tag{3.1}
\]\[
z' = -y
\tag{3.2}
\]
計算には、4次ルンゲ・クッタ法を使うことにする。\[
y_{n+1} = y_n + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4) \\
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 + \frac{1}{2}h, y_n + \frac{1}{2} k_2 h \right) \\
k_4 = f \left( x_n + h, y_n + k_3 h \right)
\tag{4}
\]
これを2関数の連立の場合\[
y' = f(x, y, z) \\
z' = g(x, y, z)
\tag{5}
\]に適用すると、\[
y_{n+1} = y_n + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4) \\
z_{n+1} = z_n + \frac{h}{6} (l_1 + 2l_2 + 2l_3 + l_4)
\tag{6}
\]となり、各係数は、\[
k_1 = f(x_n, y_n, z_n) \\
l_1 = g(x_n, y_n, z_n) \\
k_2 = f \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_1 h, z_n + \frac{1}{2} l_1 h \right) \\
l_2 = g \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_1 h, z_n + \frac{1}{2} l_1 h \right) \\
k_3 = f \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_2 h, z_n + \frac{1}{2} l_2 h \right) \\
l_3 = g \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_2 h, z_n + \frac{1}{2} l_2 h \right) \\
k_4 = f \left( x_n + h, y_n + k_3 h, z_n + l_3 h \right) \\
l_4 = g \left( x_n + h, y_n + k_3 h, z_n + l_3 h \right)
\tag{7}
\]
さて、このアルゴリズムでプログラミングしていくのですが、
以前にプログラミングの記事で書いていた
自作ウェーブ型クラスの DWave に
微分方程式求解機能をインプリメントしたら便利だろう
と思うのです。

たとえば、

double f(double x, double y);

DWave d1;
・・・ //d1 の初期設定

d1. RungeKutta4(1, f) //初期値 1, 関数 f

みたいな感じで呼び出せると便利!

と思ったのですが、考えてみたら、連立になっているので、
2値のベクトルとして扱わなければならないんですよね。

というわけで、

1.DWaeve をベクトルが扱えるように改良するか、
2.一般の連立ではなく、一関数の高階微分方程式だけに特化するか
3.DWave に実装するのはあきらめて、別でやるか。

1 は、すごく大変なので、すぐにやる気がしない^^;
2 もいいですが、とりあえず、DWave は置いておいて、3でやろうかと思っています。
数値計算>微分方程式 | コメント(0) | 2014/12/12 12:39

4次ルンゲ・クッタ法 (3) 方法の再検討

さすがに前記事の方針ではいかんともしがたいので、
方法を再検討しています。

いろいろネットで検索してみたところ、文献 [1] によると、
$f(x,y)$ の展開を y に関する展開(つまり $f_y$ や $f_{yy}$ など)だけにする
ことが可能なようです。
これができると、だいぶ楽になりますね。

ただ、その論理がどうしても理解できないで、難航しています・・・(汗)

あとは、いわゆる「4次のルンゲ・クッタ」の式はよく知られているので、
一般的な条件を導出するのはあきらめて、
その結果になるように、明らかにゼロの係数をあらかじめ消しておくというのも一案。
有名なルンゲ・クッタの式のみを導出したいのであれば、
それもありですね。

文献 [2] では、そのあらかじめ条件を制限しておく方法を取っているようです。

それにしても、この4次公式の導出は相当煩雑になるので、
皆さん、いろいろ苦労されているんだなあ・・・ということが分かりました(笑)

[1] の序文で、著者の方が

RK 法は非常に頻繁に使われる手法でありながら、
その証明の繁雑さ故に天下り的に結果を用いている場合が多く、
私もそれにならっていた。
従って講義をする立場としても内心忸怩たるものがあった。
(太字は、dyne による)

と書かれていて、気持ちわかるわーという感じです(笑)
証明をフォローしていないものを使いたくはないんですよね・・・

まだ、あきらめたわけではないですが、気長にやります^^;


参考文献
[1] 早川尚男 「Runge-Kutta法についてのノート」
http://ocw.kyoto-u.ac.jp/ja/01-faculty-of-integrated-human-studies-jp/introduction-to-simulation/pdf/rk.pdf
[2] 倭算数理研究所 「Runge-Kutta 法を導く (6) : 4次の場合」
 http://wasan.hatenablog.com/entry/2013/11/02/201200
数値計算>微分方程式 | コメント(0) | 2014/11/09 23:39

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

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

k4 を3次までテーラー展開。
\[
k_4 = f + h\left[
a_4 f_x + (b_{41} f + b_{42} k_2 + b_{43} k_3) f_y
\right] \\
+ \frac{h^2}{2} \left[
a_4^2 f_{xx}
+ 2 a_4 (b_{41} f + b_{42} k_2 + b_{43} k_3) f_{xy} \\
+ (b_{41} f + b_{42} k_2 + b_{43} k_3)^2 f_{yy}
\right] \\
+ \frac{h^3}{6} \left[
a_4^3 f_{xxx}
+ 3 a_4^2 (b_{41} f + b_{42} k_2 + b_{43} k_3) f_{xxy} \\
+ 3 a_4 (b_{41} f + b_{42} k_2 + b_{43} k_3)^2 f_{xyy} \\
+ (b_{41} f + b_{42} k_2 + b_{43} k_3)^3 f_{yyy}
\right] \\
+ \mathcal{O}(h^4)
\tag{1}
\]
この式に、前記事で求めた 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{2}
\]
および k3 の展開式
\[
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{3}
\]
を代入して、h の次数ごとに分ける。

・・・ということが果たして、人間技でできるのでしょうか?(笑)



参考文献
[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/31 12:05

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

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"
ジャンル:[学問・文化・芸術]  テーマ:[数学
数値計算>微分方程式 | コメント(0) | 2014/10/28 20:17

2次ルンゲ・クッタ法(ホイン法・中点法)

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

導きたい最終的な式を\[
y_{n+1} = y_n + h(w_1 k_1 + w_2 k_2)
\tag{1}
\]とおく。ここで、
\[
k_1 = f(x_n, y_n)
\tag{2.1}
\]\[
k_2 = f(x_n + \alpha h, y_n + \beta k_1 h)
\tag{2.2}
\]とし、$w_1$、$w_2$、$\alpha$、$\beta$ は定数とする。

$k_2$ をテーラー展開。\[
k_2 = f + \alpha h f_x + \beta k_1 h f_y + \mathcal{O}(h^2)
\tag{3}
\]以下、単に、$f$ や $f_x$ などと書いてあるのはすべて、$(x_n, y_n)$ での値を表すものとする。

(2.1)と(3)を(1)に代入。\[
y_{n+1} = y_n + h(w_1 + w_2)f + h^2 w_2 (\alpha f_x + \beta f_y f) + \mathcal{O}(h^3)
\tag{4}
\]
一方、真の解 $y(x)$ のテーラー展開は、\[
y_{n+1} = y_n + h y'(x_n) + \frac{h^2}{2} y''(x_n) + \mathcal{O}(h^3)
\tag{5}
\]ここで、元の微分方程式から\[
y'(x_n) = f
\]\[
y''(x_n) = \frac{df}{dx} = f_x + f_y y' = f_x + f_y f
\]であるから、(5) は、\[
y_{n+1} = y_n + h f + \frac{h^2}{2} (f_x + f_y f) + \mathcal{O}(h^3)
\tag{6}
\]となる。 (4) と (6) を比較して、h の2次の範囲で両者が一致するためには、\[
w_1 + w_2 = 1
\tag{7.1}
\]\[
w_2 \alpha = \frac{1}{2}
\tag{7.2}
\]\[
w_2 \beta = \frac{1}{2}
\tag{7.3}
\]であればよい。
4つの未知数に対して、3つしか条件がないから、解は無数にある。

一例として、\[
w_1 = 0, \ w_2 = 1, \ \alpha = \beta = 1/2
\]とすると、中点法の式が得られる。
\[
k_1 = f(x_n, y_n) \\
k_2 = f(x_n + h/2, y_n + hk_1/2 ) \\
y_{n+1} = y_n + h k_2
\tag{8}
\]この方式では、中点における導関数の値を評価して、
それによって次の値を求めていることになる。

また、\[
w_1 = w_2 = 1/2, \ \alpha = \beta = 1
\]とすると、ホイン法の式が得られる。
\[
k_1 = f(x_n, y_n) \\
k_2 = f(x_n + h, y_n + k_1 h) \\
y_{n+1} = y_n + \frac{h}{2} ( k_1 + k_2 )
\tag{9}
\]この方式では、積分表示\[
y_{n+1} = y_n + \int_{x_n}^{x_{n+1}} f(x,y(x)) dx \tag{10}
\]において、積分を台形近似\[
\int_{x_n}^{x_{n+1}} f(x,y(x)) dx \simeq \frac{h}{2} \{ f(x_n,y_n) + f(x_{n+1}, y_{n+1}) \}
\]により求めていることになる。

この調子で、3次、4次と進んでいくわけですが、
このままいくと、計算がどんどん複雑になっていって、
とてもできる自信がありませんね・・・^^;
途中で挫折するかもです・・・


参考文献
[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/27 22:39

1次ルンゲ・クッタ法(オイラー法)

微分方程式の数値解法について、勉強しています。

いつも、なにげなく、4次ルンゲ・クッタ公式を使って解いていますが、
この公式ってどうやって出てくるんだろうと改めて考えてみたくなりました。

まずは1次ルンゲ・クッタ法(オイラー法)から。

今後は、解きたい微分方程式を\[
y' = f(x,y)
\tag{1}
\]と仮定する。
x 軸を差分化して、$x_{n+1}$ における値 $y_{n+1} = y(x_{n+1})$ を求めるのに、
一個前の値 $(x_n, y_n)$ のみから求めることができる公式をルンゲ・クッタ型公式と呼ぶようです。

テーラー展開を用いて、\[
y_{n+1} = y_n + y'(x_n) h + \frac{y''(\xi)}{2}h^2
\tag{2}
\]ただし、$h = x_{n+1} - x_n$ で、$x_n < \xi < x_{n+1}$。

h の2次の項を落として、(1) を用いると、
1次ルンゲ・クッタ法(オイラー法)の式\[
y_{n+1} = y_n + f(x_n,y_n) h
\tag{3}
\]が得られる。

$x_n$ における導関数の値で直線近似していることに相当する。

また、(1) から直接、\[
y_{n+1} = y_n + \int_{x_n}^{x_{n+1}} f(x,y(x)) dx
\tag{4}
\]と考えることもできるので、
オイラー法は上の積分を矩形近似していることに対応している。


参考文献
[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/27 19:55
 | HOME | 

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