Octave で常微分方程式を数値計算してみましょう。最初に次のような非常に簡単な一階常微分方程式を解いてみます。
dX/dt = -X ..... (1)
X は X を微分して得られる導関数 dX/dt が X に比例(符号は負)する関数です。このような微分方程式は、放射性物質の自己崩解などでみられます。X の減少する速度が X の量に比例しているからです。これを解くために (1) 式を変型して変数分離すると次のようになります。
(1/X)dX = -1dt ..... (2)
さらに両辺を不定積分すると
log(X) = -t + C ..... (3)
C は任意の積分定数です。従って次の一般解が得られます。
X = C * exp(-t) ..... (4)
このように微分方程式の一般解は任意定数をもつので、数値積分の場合には C が一意に決まるように X の初期値 X0 = C * exp(0) が分かっている必要があります。したがって、Octave で数値積分をするときは 微分方程式の定義、tの範囲、Xの初期値 X0 の三つのデータが必要です。
Octave で常微分方程式を解くための関数は lsode 関数です。lsode は Hindmarsh の ODE solver LSODE が元になっています。lsode の書式は次のようになります。
lsode (FCN, X0, T, T_CRIT)
lsode は T に対する X の値を行列型式で出力します。FCN は Octave 上で定義された微分方程式の関数、X0 は X の初期値、T は独立変数 t の数列を納めたベクトルです。FCN は次のような型式で定義されている必要があります。
XDOT = f (X, T)
それでは、いよいよ lsode を使って式 (1) の微分方程式を解いてみましょう。つぎのように関数 xdot = f (x, t) を定義し、t の数列を作り、xの初期値 x0 = 1 で微分方程式を解きます。
octave:> function xdot = f (x, t)
> xdot = -x; .... これは常微分方程式 dx/dt = - x の式そのままです。
> endfunction
octave:> t = linspace (0, 5, 100);
octave:> x0 = 1;
octave:> y = lsode ("f", x0, t);
octave:> plot (y)
予想通り指数関数になっているようです。
微分方程式の数値計算法のやり方は概ね次のようになります。最初に x の初期値 x0 と t0 = 0 から計算を始めます。xdot = f(x, t) という関数からは x = x0, t = t0 の時の dx/dt を計算することができます。すると t = t1 のときの x の値 x1 は x1 = x0 + (dx/dt|t=t0) * Δt で計算できます。この新しい x1 と t1 を使って次の dx/dt|t=t1 を計算することによって x2 を得ることができます。これをくり返すことによって x0, x1, x2, ...., xn を計算することができるのです。もちろんこのオイラー法そのままでは精度が悪いので lsode では精度を上げる工夫をしてありますが、基本的な考え方はおなじです。
次に述べる高階の常微分方程式も変数をベクトルにすることによって高階常微分方程式を一階常微分方程式に変換して計算をおこなうので基本的には一階常微分方程式と同じアルゴリズムで数値計算が行われます。
今度は次のような二階常微分方程式を解いてみましょう
x'' + x' + 2*x = 1 ..... (5)
lsode は dX/dt = f(x, t) という型の常微分方程式しか解くことができませんから、式 (5) はそのままでは解くことができません。そこで x1 = x, x2 = dx/dt となるように変数を置き換えます。すると式 (5) は次のように変型できます。
dx1/dt = x2 ..... (6)
dx2/dt = -x2 - 2*x1 + 1 .... (7)
ここで X = [x1; x2] とおくと、(6)(7) の連立方程式は dX/dt = f(X, t) とベクトル変数関数の一階常微分方程式になっているのが分かります。それでは式 (5)を初期値 X0 = [0; 0] の条件の元で解いて見ましょう。(この場合 t = 0 のとき x(0) = 0, x'(0) = 0 になります。)
octave:> function xdot = f (x, t)
> xdot(1) = x(2);
> xdot(2) = -x(2) - 2*x(1) + 1;
> endfunction
octave:> t = linspace (0, 20, 100);
octave:> x0 = [0; 0];
octave:> y = lsode ("f", x0, t);
octave:> plot (y)
plot で描かれるグラフの line 1 は x(t) を示し、line 2 は x'(t) を示しています。
二階常微分方程式の応用で等加速度運動を調べてみましょう。等加速度運動の運動方程式は
x'' = a
です。これを lsode で解くと次のようになります。
octave:> function xdot = f (x, t)
> xdot(1) = x(2);
> xdot(2) = 1;
> endfunction
octave:> t = linspace (0, 5, 100);
octave:> x0 = [0; 0];
octave:> y = lsode ("f", x0, t);
octave:> plot (y)
plot の結果から速度は直線的に、距離は時間の二乗に比例して増加していることが分かります。
壁にバネ定数 k のバネで固定された質量 m の物体の水平方向の振動を考えてみましょう。物体 m には粘性係数 c の粘性減衰がかかっているとします。物体 m の加速度が x'' とすると物体の運動エネルギーは m * x'' です。また m には速度に比例する減衰力がかかりますがこれは c * x' です。また m にはバネの復元力 k * x がかかります。これらが全体として m にかかる外力 f(x) とつりあっているのでこの力学モデルの運動方程式は
m * x'' + c * x' + k * x = f(t)
となります。これは二階常微分方程式なので一階の連立常微分方程式になおします。
x1' = x2
x2' = -(c/m)*x2 - (k/m)*x1 + f(t)
ここで x1 = 1, x1' = 0 を初期値として、外力 f(t) が全くかからない場合の m の動きをみるスクリプト vibration.m を次のように作成します。
$ cat vibration.m
global m c k;
m = input ("m = ");
c = input ("c = ");
k = input ("k = ");
function xdot = f (x, t)
global m c k;
xdot(1) = x(2);
xdot(2) = -(c/m)*x(2) - (k/m)*x(1);
endfunction
x0 = [1; 0];
t = linspace(0, 25, 100);
y = lsode("f", x0, t);
plot(t, y(:,1))
このスクリプトを使って、m c k の値をいろいろと変えたときの物体 m の運動を観察します。まず m = 1, c = 0.4, k = 1 で試してみます。これは減衰振動になります。
$ octave GNU Octave, version 2.0.16 (i386-vine-linux-gnu). Copyright (C) 1996, 1997, 1998, 1999, 2000 John W. Eaton. This is free software with ABSOLUTELY NO WARRANTY. For details, type `warranty'. octave:1> source "vibration.m" m = 1 c = 0.4 k = 1
粘性係数を 0 にすると単振動になります。
octave:2> source "vibration.m" m = 1 c = 0 k = 1
質量を減らすと振動数が増加します。
octave:3> source "vibration.m" m = 0.5 c = 0 k = 1
粘性係数を上げると無周期減衰運動になります。
octave:4> source "vibration.m" m = 1 c = 2 k = 1