情報基礎「Pythonプログラミング」（ステップ６：常微分方程式）

Pythonプログラミング（ステップ６・常微分方程式）

このページの目標

常微分方程式の数値解法が反復計算（総和のパターン）に還元できることを理解すること。
簡単な反復解法をコーディングできること。

１階の常微分方程式の数値解法

一階の常微分方程式（ordinary differential equation) $$ \frac{d x(t)}{d t} = f(t, x(t)) \tag{1} $$ を考えよう。$x$は「座標」、$t$は「時刻」を表すような変数であると想定して説明することにする。

もしも、右辺が$t$のみの関数であれば ($f(t)$ ならば)、この解は、初期条件 $x(t_0) = x_0$ のもとで、 $$ x(t) = x_0 + \int_{t_0}^t f(t) dt \tag{2} $$ と書けるから、$f(t)$の積分が計算できれば、微分方程式は解かれたことになる。もし$f(t)$の積分が解析的に求まらない場合でも、台形法やシンプソン法などを用いて数値的に積分を評価することで、任意の $t$ での $x(t)$ が（近似的ながら）得られる。

ところが、式(1)の右辺に $x(t)$ が含まれると、式(2)のように直接的に解を表現することはできない。そこで、初期値から出発して、次々と「次の時刻」の $x$ を求める（予測する）、というアプローチが採られる。中でも、最もシンプルなのが、以下に述べるオイラー法(Euler method)である。

そもそも、微分は差分の極限として定義されたことを思い出すと、以下の式で $\Delta t$ が十分小さければ、微分は差分でよく近似できる。すなわち、 $$ \frac{d x(t)}{dt} \approx \frac{x(t+\Delta t) - x(t)}{\Delta t} $$ である。そこで、式(1) の左辺を差分で置き換え、 $$ \frac{x(t+\Delta t) - x(t)}{\Delta t} = f(t, x(t)) $$ と置いてみる。この式を少し変形すると $$ x(t + \Delta t) = x(t) + f(t, x(t)) \Delta t \tag{3} $$ となる。ここで、式(3)の右辺には「時刻」$t$ での値のみが登場している。つまり、時刻 $t$ の状態を使って、次の時刻 $t + \Delta$ での $x$ の値が与えられる。このように、既知の値のみを使って、次のステップの値を計算する方式を陽解法(explicit method)と呼ぶ。

陰解法

一方、微分を $$ \frac{d x(t)}{dt} \approx \frac{x(t) - x(t - \Delta t)}{\Delta t} $$ のような差分（後退差分）で近似しても良いので、微分方程式は $$ \frac{x(t) - x(t - \Delta t)}{\Delta t} = f(t, x(t)) $$ のように書くこともできる。すると（時間を$+\Delta t$だけずらすと）、式(3)の代わりに、 $$ x(t + \Delta t) = x(t) + f\left(t + \Delta t, x(t+\Delta t) \right) \Delta t \tag{4} $$ としても良いはずである。この式は右辺にも $x(t+\Delta t)$ が含まれているから、「次の時刻の値が、次の時刻の値にも依存する」という、自己言及的な構造になっている。このように、未知の値を使って未知の値を計算する方式は、陰解法(implicit method)と呼ばれる。陰解法で微分方程式を解くには、各ステップで、$x(t +\Delta t)$ を未知数とした方程式を解く手間が生じるが、安定性の点では一般に有利である。

以上を踏まえ、初期条件$x(t_0)=x_0$から出発し $t \le T$ まで、陽解法のオイラー法で常微分方程式（１変数）を解くためのアルゴリズムをまとめると、次のようになる：

$t \leftarrow t_0$
$x \leftarrow x_0$
$\Delta t$ に「十分に小さな」数値をセットする
$t \lt T$ の間、以下を反復する：
　$x \leftarrow x + f(t,x) \Delta t$
　$t \leftarrow t + \Delta t$
反復ここまで

このアルゴリズムは、基本的には総和計算と全く同じパターンであることに気づくはずだ。具体的に、$\frac{dx}{dt} = -x$ をオイラー法（陽解法）で解くコードの例を以下に示す。初期条件は $x(0)=1.0$ とした。

exp-decay.py

この微分方程式の数学的な解は $x(t)=\exp(-t)$ であるので、数値解と比較できるようコードを修正してみると良い。

# coding: utf-8

dt=0.01
T=2.0

t=0.0
x=1.0
while t<T:
    print("t=",t," x=",x)
    f=-x
    x=x+f*dt
    t=t+dt

　練習：微分方程式

オイラー法（陽解法）を用いて、二つのパラメータ（$r, N$)を持つ微分方程式 $$ \frac{dx}{dt} = r x \left(1 - \frac{x}{N} \right) $$ の数値解を求めよ。初期条件は$x(0)=1$とせよ。また、$r=1.0, N=100$ とする。

　ヒント

この微分方程式はロジスティック方程式として知られ、生物の個体数や人口の変化を考える際の基本的なモデルである。パラメータ$r$は個体の自発的な増殖率を、$N$は「上限」を与える。 $r$と$N$を変えながら、変化の様子をシミュレーションしてみなさい。

2階の常微分方程式の数値解法

物理学などで登場する微分方程式の多くは二階微分を含む。例えば、長さ $L$ の軽い棒と錘で構成された振り子の運動は、振れ角を$\theta$とすると、微分方程式 $$ \frac{d^2 \theta}{d t^2} = - \frac{g}{L} \sin(\theta) $$ で記述される。ここで $g$ は重力加速度である。

このような二階の微分方程式は、全く等価な、２つの連立した一階の微分方程式に書き直すことができる。振り子の例では、$x(t) = \theta(t)$, $y = \frac{d \theta(t)}{dt}$ と改めて変数を置き直すと、 $$ \begin{align} \frac{dx}{dt} & = y \\ \frac{dy}{dt} & = - \frac{g}{L} \sin(x) \\ \end{align} $$ となる。

それぞれの方程式にオイラー法を適用すると、以下のようなコードを書くことができるだろう。

# coding: utf-8
import math

g=9.8
L=1.0
dt=0.01
T=3.0

t=0.0
x=0.2
y=0.0
while t<T:
    print("t=",t," x=",x," y=",y)
    x1=x+y*dt
    y1=y-g/L*math.sin(x)*dt
    x=x1
    y=y1
    t=t+dt

ここで、$g=9.8 \ \textrm{m}/\textrm{s}^2$, $L = 1\ \textrm{m}$ とし、初期状態の触れ角は $0.2\ \textrm{rad}$, 初期角速度は $0 \ \textrm{rad}/\textrm{s}$ とおいた。

上記のコードで、変数 x1, y1は、次の時間ステップの値を記憶しておくために使っており、両方の値が求まった時点で、改めて x, yを更新するようになっている。微分方程式を数値的に扱う際には、変数に「いつの時点」の値が保持されているのかを注意深く考える必要がある。

エネルギー保存を破らないように数値計算するためには、オイラー法よりも高精度の公式（例えばこちら）を用いたり、保存系の計算に適した方法（シンプレクティック解法）を選ぶ必要がある。

このコードを実行すると、振動の様子が再現されるかに見えるが、計算する時間を長くすると(変数Tの値を大きくすると）、振幅と角速度が徐々に増大することに気づく。すなわち、本来ならば保存されるべき系の力学エネルギーが増加するという、非現実的な挙動を示す。これは、振り子の系にオイラー法を適用すると、反復の過程で誤差が（一定のバイアスを持って）累積してしまうためで、一見正しそうなデータが得られた場合でも、その結果をよく吟味する必要がある。

　練習：摩擦のある振り子

摩擦があるような振り子の系 $$ \frac{d^2 \theta}{d t^2} = - \frac{g}{L} \sin(\theta) - \gamma \frac{d\theta}{dt} $$ をオイラー法で計算してみなさい。ここで $\gamma \ [1/\textrm{s}]$ は回転に伴う摩擦の係数（ここでは、振り子の錘が速度に比例した抵抗力を受けると仮定。質量 $m$ の錘が速度 $V$ で運動する際の抵抗力を $\eta V$ とすると、$\gamma = \eta/m$）。

$g=9.8 \ \textrm{m}/\textrm{s}^2$, $L = 1\ \textrm{m}$ とし、初期状態の触れ角は $0.2\ \textrm{rad}$, 初期角速度は $0 \ \textrm{rad}/\textrm{s}$ とし、 $\gamma$ の値を変えて摩擦の効果をシミュレーションしてみること。