連立微分方程式
前回は1階の微分方程式を解く方法について解説しました。しかし実際に出会う問題は、多くの場合2階の微分方程式になっています（例えば運動方程式）。そこで今回は、高階の微分方程式の解法について考えましょう。よく知られているように、高階の微分方程式は連立の微分方程式に帰着できます。そこではじめに、連立微分方程式の解法を解説します。
次のような連立微分方程式を考えます。
\left.
\begin{array}{l}
\frac{dy_1(x)}{dx} &=& f_1(x, y_1,\cdots , y_n) \\
\frac{dy_2(x)}{dx} &=& f_2(x, y_1,\cdots , y_n) \\
\cdots \\
\frac{dy_n(x)}{dx} &=& f_n(x, y_1,\cdots , y_n) \\
\end{array}
\right\}\\[35pt]
y_1(x_0) = y_1^{(0)},\ y_2(x_0) = y_2^{(0)}, \cdots , y_n(x_0) = y_n^{(0)}
これをベクトルを用いて表現すると次のようになります。
\left.
\begin{array}{l}
\frac{d{\bf y}(x)}{dx} &=& {\bf f}(x, {\bf y}(x)) \\
{\bf y}(x_0) &=& {\bf y^{(0)}}
\end{array}
\right\}
ここで
{\bf y}(x) = \left(
\begin{array}{l}
y_1(x) \\
y_2(x) \\
\vdots \\
y_n(x)
\end{array}
\right), \ \ \ \ \ {\bf f}(x, {\bf y}(x)) = \left(
\begin{array}{l}
f_1(x,{\bf y}(x)) \\
f_2(x,{\bf y}(x)) \\
\vdots \\
f_n(x,{\bf y}(x))
\end{array}
\right), \ \ \ \ \ {\bf y^{(0)}} = \left(
\begin{array}{l}
y_1^{(0)} \\
y_2^{(0)} \\
\vdots \\
y_n^{(0)}
\end{array}
\right)
です。この連立常微分方程式を解く方法は、基本的に連立でない場合と同様です。
それでは次の連立常微分方程式を4次のルンゲ・クッタ法で解いてみましょう。
\left.
\begin{array}{l}
\frac{dy_1}{dx} &=& y_2 \\
\frac{dy_2}{dx} &=& -1
\end{array}
\right\}
ただし、
$$y_1(0) = 0, \ \ \ \ \ y_2(0) = 1$$
とします。この連立微分方程式の厳密解が
\left.
\begin{array}{l}
y_1 &=& x - \frac{1}{2}x^2 \\
y_2 &=& 1 - x
\end{array}
\right\}
となることは簡単に確かめられるでしょう。
#include <stdio.h>
#include <stdlib.h>
double f(int i, double x, double *y);
void runge4(double x, double *y, int n, double h, double funcpt(int i, double x, double *y));
int main(void)
{
double x, y[2], h;
int n, imax, i;
x = 0.0;
y[0] = 0.0;
y[1] = 1.0;
n = 2;
h = 0.1;
imax = 2.0 / h + 1;
printf("x=%12.4e y1=%12.4e%12.4e y2=%12.4e%12.4e\n", x, y [0], x - 0.5 * x * x, y[1], 1.0 - x);
for (i = 1; i <= imax; i++) {
runge4(x, y, n, h, f);
x = x + h;
printf("x=%12.4e y1=%12.4e%12.4e y2=%12.4e%12.4e\n", x, y [0], x - 0.5 * x * x, y[1], 1.0 - x);
}
return 0;
}
double f(int i, double x, double *y)
{
if (i == 0) return y[1];
else if (i == 1) return -1.0;
else exit(1);
}
/* 4次のルンゲ・クッタ法により連立微分方程式を解くための関数 */
void runge4(double x, double *y, int n, double h, double funcpt(int i, double x, double *y))
{
double *k1pt, *k2pt, *k3pt, *k4pt, *work;
int i;
k1pt = (double *) malloc(n*8);
k2pt = (double *) malloc(n*8);
k3pt = (double *) malloc(n*8);
k4pt = (double *) malloc(n*8);
work = (double *) malloc(n*8);
for (i = 0; i < n; i++) k1pt[i] = h * funcpt(i, x, y);
for (i = 0; i < n; i++) work[i] = y[i] + k1pt[i] / 2.0;
for (i = 0; i < n; i++) k2pt[i] = h * funcpt(i, x + h / 2.0, work);
for (i = 0; i < n; i++) work[i] = y[i] + k2pt[i] / 2.0;
for (i = 0; i < n; i++) k3pt[i] = h * funcpt(i, x + h / 2.0, work);
for (i = 0; i < n; i++) work[i] = y[i] + k3pt[i];
for (i = 0; i < n; i++) k4pt[i] = h * funcpt(i, x + h, work);
for (i = 0; i < n; i++) y[i] = y[i] + (k1pt[i] + 2.0 * k2pt[i] + 2.0 * k3pt[i] + k4pt[i]) / 6.0;
free(k1pt);
free(k2pt);
free(k3pt);
free(k4pt);
free(work);
}
C 言語では、プログラムの実行中に新たな記憶領域を確保できます（ malloc ）。したがって、作業用の配列や変数を必要の応じて関数内で確保できるため、関数を呼ぶ側（今の場合はメイン関数）で配列を宣言して引数で関数に渡すといった面倒なことをする必要はありません。なお malloc 関数を利用する場合、<stdlib.h> のヘッダファイルをインクルードする必要があります。
課題
次の連立常微分方程式を解きなさい。
\left.
\begin{array}{l}
\frac{dy_1}{dx} &=& y_2 \\
\frac{dy_2}{dx} &=& -y_1
\end{array}
\right\}
ただし、
$$y_1(0) = 0, \ \ \ \ \ y_2(0) = 1$$
とします。また厳密解を求め、それと比較しなさい。