FC2ブログ

二階微分方程式の数値求解 (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
コメント

管理者のみに表示