Maxima で常微分方程式

常微分方程式の解法 ode2()

Maxima は 1 階および 2 階の常微分方程式を ode2 関数を使って解くことができます。

Maxima を使って常微分方程式を解く例については，理工系の数学C でやりましたので，詳細は省略します。以下を参照してください。

• 理工系の数学C
  • Maxima-Jupyter で常微分方程式

練習 2-1

重力場中の投げ上げ運動

1. 高さ $0$ から初速度 $v_0$ で鉛直上方に投げ上げた物体の時刻 $t$ における高さ $y$ を求めなさい。運動方程式は，
   $$ \frac{d^2 y}{dt^2} = -g$$
2. 速度に比例する空気抵抗がある場合，運動方程式は以下のようになる。
   このとき，高さ $0$ から初速度 $v_0$ で鉛直上方に投げ上げた物体の時刻 $t$ における高さ $y$ を
   求めなさい。
   $$ \frac{d^2 y}{dt^2} = -g – \beta \frac{dy}{dt}$$
3. 前問 2. で求めた解が $\beta \rightarrow 0$ のとき，前問 1. の答えに一致することを示しなさい。

ode2() で解いた微分方程式の解に含まれる $\%k_1, \%k_2$ は積分定数です。（2階微分方程式の解には積分定数が2個。）

$t=0$ のとき $\displaystyle y = 0, \ \frac{dy}{dt} = v_0$ という初期条件を課して積分定数を決める例です。

$t=0$ のとき $\displaystyle y = 0, \ \frac{dy}{dt} = v_0$ という初期条件を課して積分定数を決める例です。

空気抵抗がある場合の運動方程式の解は分母に $\beta$ があり，単純に $\beta = 0$ にすると大変なことになってしまいそうです。

そこで，以下のように極限をとってみます。

$$\lim_{\beta \rightarrow 0} y$$

常微分方程式の数値解法 rk()

常微分方程式が解析的に解けないときは，数値的解法によって近似値を求めることができます。Maxima にはルンゲ・クッタ法によって常微分方程式を数値的に解く関数 rk() があらかじめ組み込まれています。使い方を示します。

1階常微分方程式

1階の常微分方程式
$$ \frac{dx}{dt} = f(x, t)$$
をルンゲ・クッタ法で解く例です。

初期条件を $t=t_0$ で $x(t_0) = x_0$ とし，$t=t_0$ から $t = t_1$ までを $N$ 分割し，刻み幅 $\displaystyle h = \frac{t_1 – t_0}{N}$ で $x$ について解く場合，以下のようにします。

rk(f(x,t), x, x0, [t, t0, t1, h])

以下では，簡単なロジスティック方程式
$$ \frac{dx}{dt} = (1-x)\,x$$
を，初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 0.1$，$t = t_0 = 0$ から $t_1 = 10$ までを $ N = 100$ 分割して解きます。

Maxima の rk() はシンプルな4次のルンゲ・クッタ法ですから，きざみ幅 h の取り方が誤差に直結します。数値計算の誤差については，h を変えて計算してみて精度を確かめる必要があります。

ルンゲ・クッタ法による数値解は，リスト rkdat に代入しています。リスト rkdat の「長さ」（項数）は $N+1$ 個です。

リスト rkdat の最初の項 rkdat[1] を表示させると，$[t_0, x_0]$，最後の項 rkdat[length(rkdat)] の値は $[t_1, x(t_1)]$ です。

刻み幅を変更すると数値解の精度が変わります。念のため，$N$ を変更して数値解を確認します。

結果を見比べた限りでは，数値解は $N=400$ とした場合に
$$ x(t_1) = 0.9995915675…$$ あたりまでの精度がありそうです。

ルンゲ・クッタ法による数値解は，リスト rkdat に代入されます。
これをグラフにしてみます。

2階常微分方程式

次のような2階常微分方程式をルンゲ・クッタ法で解く。

$$ \frac{d^2 x}{dt^2 } = f\left(x, \frac{dx}{dt}, t\right)$$

この場合には， $\displaystyle v \equiv \frac{dx}{dt}$ とおけば，次のような連立1階常微分方程式の形に帰着できる。
\begin{eqnarray}
\frac{dx}{dt} &=& F_1(x, v, t) = v \\
\frac{dv}{dt} &=& F_2(x, v, t) = f(x, v, t)
\end{eqnarray}

初期条件を $t=t_0$ で $x(t_0) = x_0, v(t_0) = v_0$ とし，$t=t_0$ から $t = t_1$ までを $N$ 分割し，刻み幅 $\displaystyle h = \frac{t_1 – t_0}{N}$ で $x, v$ について解く場合，以下のようにします。

rk([F1(x,v,t), F2(x,v,t)], [x, v], [x0, v0], [t, t0, t1, h])

以下では例として，簡単な減衰＋強制振動の方程式

$$ \frac{d^2 x}{dt^2} = -x – a \frac{dx}{dt} + b \cos(t) $$

を，

• $a , b $ にいろいろな値を入れて
• 初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 3, v_0 = v(t_0) = 0$，
• $t = t_0 = 0$ から $t_1 = 20$ までを $ N $ 分割

して解きます。

上記の結果から，だいたい $h = 0.01$ の刻み幅で小数点以下8桁程度の精度があると推測されます。以下では，この値で計算を続けます。

$a = 0, b = 0$ の場合の数値解を dataa00b00，$a = 0.5, b = 0.2$ の場合の数値解を data05b02 に代入します。

以下の例からわかるように，これらのリストは1列目が $t$, 2列目が $x$, 3列目が $v$ の値です。

1列目のデータを横軸，2列目のデータを縦軸にとってグラフを描く例です。

makelist([c[1],c[2]], c, data00b00) がキモです。

2つのデータのグラフを重ねて表示する例です。凡例やグリッドも設定します。

練習 2-2

$ a = 0.5, b = 0$ の減衰振動の場合の数値解を同様にしてもとめ，
上記の $a = 0.5, b = 0.2$ の減衰＋強制振動の解と比較せよ。

参考：減衰＋強制振動の解析解

実は Maxima では，以下のような減衰＋強制振動の微分方程式を ode2() によって解析的に解くことができます。

以下では，
$$ \frac{d^2 x}{dt^2} = -x – a \frac{dx}{dt} + b \cos(t) $$
を，$a = 0.5, b = 0.2$ の場合に解析的に解く例を示します。

解くべき微分方程式の記述例です。微分 diff() の前に ' をつけて書きます。

上記の解は，2つの積分定数 $\%k_1, \%k_2$ を含む一般解です。

2階常微分方程式の初期条件は ic2() 関数で与えます。以下では，初期条件として $t = 0$ で $x = 3, \frac{dx}{dt} = 0$ を設定しています。

この解を関数 x(t) として定義します。

ルンゲ・クッタ法 rk() で求めた数値解と一緒に図示します。数値解は2001個と多いので，間引いたリスト rkdat を作り，描画スタイルを points にしてグラフにします。

念のために，$t=t_1$ での数値解と解析解を比べてみます。 （小数点以下9桁程度の精度があることがわかるでしょうか。）