常微分方程式
次のような1階の常微分方程式の解法を考えます。
$$\frac{dy(x)}{dx} = f(x, y(x)) \tag{1}$$
ここで初期条件を $y(x_0) = y_0$ とします。初期条件が与えられたときに微分方程式を解く問題を、初期値問題とよびます。初期値問題の解法にはいろいろありますが、ここではオイラー法とルンゲ・クッタ法を取り上げます。
オイラー法
まず式 (1) において与えられた初期条件 $y(x_0) = y_0$ から、刻み幅 $h$ だけ離れた点 $x_1 = x_0 + h$ における $y(x)$ の値を求めましょう。$y(x_1) = y(x_0 + h)$ を $x = x_0$ のまわりでテイラー展開します。
$$y(x_1) = y(x_0) + y'(x_0)h + O(h^2)$$
ここで $y'(x_0) = f(x_0, y_0)$ なので、上式は
$$y(x_1) = y(x_0) + hf(x_0, y_0) + O(h^2)$$
となります。$h^2$ のオーダーを無視して $y(x_1)$ の近似解 $y_1$ を
$$y_1 = y_0 + hf(x_0, y_0)$$
として求めます。同様にして $x = x_2 = x_1 + h$ における近似解は
$$y_2 = y_1 + hf(x_1, y_1)$$
となります。一般に $x = x_k$ における近似解は
$$y_k = y_{k-1} + hf(x_{k-1}, y_{k-1})\ \ \ \ \ \ \ \ (k = 1,2,\cdots)$$
で与えられます。
それでは、次の微分方程式をオイラー法で解いてみましょう。
$$\frac{dy}{dx} = -y\ \ \ \ \ \ \ \ (\ y(0) = 1\ )$$
この微分方程式の厳密解はもちろん $y=e^{-x}$ です。
#include <stdio.h>
#include <math.h>
double f(double x, double y);
void euler(double x, double *y, double h, double funcpt(double x, double y));
int main(void)
{
double x, y, h;
int i, imax;
x = 0.0;
y = 1.0;
h = 0.05;
imax = 1.0 / h + 1;
printf("x=%12.4e y=%12.4e%12.4e\n", x, y, exp(-x));
for (i = 1; i <= imax; i++) {
euler(x, &y, h, f);
x = x + h;
printf("x=%12.4e y=%12.4e%12.4e\n", x, y, exp(-x));
}
return 0;
}
double f(double x, double y)
{
return -y;
}
/* オイラー法により常微分方程式を解くための関数 */
void euler(double x, double *y, double h, double funcpt(double x, double y))
{
*y = *y + h * funcpt(x, *y);
}
2次ルンゲ・クッタ法
まず一番簡単な2次の公式を説明します。$y(x_1)$ のテイラー展開を $h^2$ まで行います。
\begin{eqnarray}
y(x_1) &=& y(x_0) + y'(x_0)h + \frac{1}{2!}y''(x_0)h^2 + O(h^3) \\
&=& y_0 + hf(x_0, y_0) + \frac{1}{2!}h^2f'(x_0, y_0) + O(h^3) \\
&=& y_0 + hf(x_0, y_0) + \frac{1}{2!}h^2[f_x(x_0, y_0) + f_y(x_0, y_0)f(x_0, y_0)] + O(h^3) \tag{2}
\end{eqnarray}
$h^3$ 以上を無視してこの公式をそのまま使うと、$f(x, y)$ の偏微分を計算しなければなりません。これを避けるために、次のような形の偏微分を含まない近似解を考えます。
$$y_1 = y_0 + ahf(x_0, y_0) + bhf(x_0 + ph, y_0 + qk_1)$$
ここで $a, b, p, q$ は未定のパラメータです。この式の右辺の $f(x_0 + ph, y_0 + qk_1)$ を再びテイラー展開すると
$$y_1 = y_0 + h(a+b)f(x_0, y_0) + h^2b[pf_x(x_0, y_0) + qf_y(x_0, y_0)f(x_0, y_0)] + O(h^3)$$
となります。この式と式 (2) とが $h^2$ の項まで一致するようにパラメータを選びます。すなわち
$$a + b = 1,\ \ \ \ \ \ \ \ bp = bq = \frac{1}{2}$$
とします。結局自由に選べるパラメータは1つだけになりました。一般的な公式は、ここでは $b$ をパラメータとして
\begin{eqnarray}
y_k &=& y_{k-1} + (1-b)k_1 + bk_2 \\
k_1 &=& hf(x_{k-1}, y_{k-1}) \\
k_2 &=& hf(x_{k-1} + \frac{h}{2b}, y_{k-1} + \frac{k_1}{2b})
\end{eqnarray}
と書けます。このように、近似解のテイラー展開が $h^2$ の項まで真の解のテイラー展開と一致するという意味で、これを2次の公式とよびます。$b$ の値は $0$ 以外ならば任意に選べますが、普通は $b = 1/2$ や $b = 1$ 等が用いられるようです。
オイラー法では、 $x = x_1$ での近似解を求めるのに $x = x_0$ における $f(x, y)$ の値しか使いませんでした。しかしここでは、$(x, y) = (x_0 + h/2b, y_0 + k_1/2b)$ での $f(x, y)$ の値も使い、両者の重み付き平均を使って近似解を求めています。このように $x = x_0$ 以外での $f(x, y)$ の値を用いて近似解をを計算する公式をルンゲ・クッタ型公式とよびます。
それでは2次の公式を使って、先ほどと同じ問題を解いてみましょう。
#include <stdio.h>
#include <math.h>
double f(double x, double y);
void runge2(double x, double *y, double h, double b, double funcpt(double x, double y));
int main(void)
{
double x, y, b, h;
int i, imax;
x = 0.0;
y = 1.0;
h = 0.05;
b = 0.5;
imax = 1.0 / h + 1;
printf("x=%12.4e y=%12.4e%12.4e\n", x, y, exp(-x));
for (i = 1; i <= imax; i++) {
runge2(x, &y, h, b, f);
x = x + h;
printf("x=%12.4e y=%12.4e%12.4e\n", x, y, exp(-x));
}
return 0;
}
double f(double x, double y)
{
return -y;
}
/* 2次のルンゲ・クッタ公式により常微分方程式を解くための関数 */
void runge2(double x, double *y, double h, double b, double funcpt(double x, double y))
{
double k1, k2;
double y1;
k1 = h * funcpt(x, *y);
y1 = *y + k1 / (2.0 * b);
k2 = h * funcpt(x + h / (2.0 * b), y1);
*y = *y + (1.0 - b) * k1 + b * k2;
}
4次ルンゲ・クッタ法
さて、オイラー法と2次のルンゲ・クッタ法のプログラムを実際に実行した人は気がついたかと思いますが、これらの解法はあまり精度がよくありません。これらの解法は単純なので、微分方程式を数値的にどのようにして解くのかという解法の原理を理解するのには適しています。しかし、実際に微分方程式を解くときにオイラー法や2次の公式が用いられることは普通はありません。そこで最後に、ルンゲ・クッタ法でよく用いられる4次の公式について紹介します。この公式は2次の公式を導出した方法と同様に、解を $h^4$ の項までテイラー展開することによって得られます。ここではその導出は省略し、結果だけを示します。
\begin{eqnarray}
y_k &=& y_{k-1} + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \\
k_1 &=& hf(x_{k-1}, y_{k-1}) \\
k_2 &=& hf(x_{k-1} + \frac{h}{2}, y_{k-1} + \frac{k_1}{2}) \\
k_3 &=& hf(x_{k-1} + \frac{h}{2}, y_{k-1} + \frac{k_2}{2}) \\
k_4 &=& hf(x_{k-1} + h, y_{k-1} + k_3)
\end{eqnarray}
余力のある人は4次の公式の導出にもトライしてみてください。ところで、4次の公式の場合には自由に選べるパラメータが2つ存在します。このパラメータの値を変えることによってルンゲ・クッタ・ギル法とよばれる公式が得られます。数値計算の本やインターネット上の解説などを参考に、この公式についても学んでおくとよいでしょう。
課題
次の微分方程式をオイラー法で解きなさい。
$$\frac{dy}{dx} = -xy + y^2\ \ \ \ \ \ \ \ (\ y(1) = 1\ )$$
また次の微分方程式をルンゲクッタ法で解きなさい。解は $x = 0$ から $x = \pi/2$ まで求めなさい。
$$\frac{dy}{dx} = \sin(x) + \cos(y)\ \ \ \ \ \ \ \ (\ y(0) = 0\ )$$
最後に4次の公式を用いたプログラムを完成させ、いろいろな微分方程式を実際に解きなさい。