Webページでは数式を表すためにLaTeX表記が使えるMathJaxを利用しています。 WebブラウザにはSafari/Chrome/Firefoxを使って下さい(IEでは表示できないようです。)

微分方程式の数値解法

$\boldsymbol{x}(t)=(x_1(t),\dots,x_n(t))$ の初期条件 $\boldsymbol{x}(t_0)=(x_1(t_0),\dots,x_n(t_0))$ のもとでの微分方程式系 \begin{align*} \frac{dx_1}{dt}&=f_1(\boldsymbol{x}, t), \quad x_1(t_0)=x_{10}\\ \vdots &\\ \frac{dx_n}{dt}&=f_n(\boldsymbol{x}, t), \quad x_n(t_0)=x_{n0} \end{align*} をベクトル表記して、次のように表す。 \[ \frac{d\boldsymbol{x}}{dt}=\boldsymbol{f}(\boldsymbol{x},t),\qquad \boldsymbol{x}(t_0)=\boldsymbol{x}_0 \] この微分方程式は各点 $\boldsymbol{x}$ および時間 $t$ 毎に接線ベクトル $\boldsymbol{f}$ を与えているとみることができ、$\boldsymbol{f}$ をベクトル場ということがある。

時間 $\Delta t$ の経過による $\boldsymbol{x}$ の変化を考えると \[ \lim_{\Delta t\rightarrow 0}\frac{\boldsymbol{x}(t+\Delta t)-\boldsymbol{x}(t)}{\Delta t}= \frac{d\boldsymbol{x}}{dt} \] であることから、有限ではあるが十分小さい $\Delta t$ に対しては \[ \frac{\boldsymbol{x}(t+\Delta t)-\boldsymbol{x}(t)}{\Delta t}\approx \boldsymbol{f}(\boldsymbol{x},t) \] が期待できる。 そこで時間刻み $\Delta t$ ごとに、$t$ をパラメータとする解曲線 $\boldsymbol{x}(t)$ 上に近いと考えられる離散的な値を逐次的に求める微分方程式の数値解法がさまざまに考案されてきた。

ここでは簡単のために、時間刻みを $\Delta t$ を一定値と考える。 \[ \boldsymbol{x}_n\equiv \boldsymbol{x}(t_0+n\Delta t),\qquad n=0,1,2,\dots \] と表記する。

Euler法

\[ \boldsymbol{x}_{n+1}= \boldsymbol{x}_n +\Delta t \boldsymbol{f}(\boldsymbol{x}_n, t_n),\quad t_n = t_0 + n\Delta t \] にしたがって点列 $\boldsymbol{x}_0, \boldsymbol{x}_1,\dots, \boldsymbol{x}_n,\dots$ を求める方法をEuler法という。

次の関数 euler(x0, T, dt, f) は、Euler法を使って、時刻 $t_0$ における$m$次元(2行目の width に次元の値が格納される)の初期点リスト x0 =[$x_{01}$,$\dots$,$x_{0m}$]から出発して、時刻 $t_0 + T$ までを指定した時間刻み dt によって $N=T/dt$ 個からなる点のリスト [$\boldsymbol{x}_0$, $\boldsymbol{x}_1$, $\dots$, $\boldsymbol{x}_N$] を返す(3行目で $\boldsymbol{x}_0$ を $\boldsymbol{x}$ とおいて、while文で更新した点を orbit.append(list(x)) によって数値軌道リスト orbit に次々と追加している)。 f はベクトル場 $\boldsymbol{f}$ を定める $m$ 個の関数リスト [$f_1$, $\dots$, $f_m$] である。

Python的技法は、7行目の fx = map(lambda a: a(t, *x), f) である(lambda a: a(t, *x) がラムダ計算による関数定義で、その引数として t, *x をセットしている。下の可変長引数)参照。 Python関数 map と lamda計算を使って、ベクトル $\boldsymbol{f}(t,\boldsymbol{x})$ をリスト fx = [$f_1(\boldsymbol{x}, t)$, $\dots$, $f_m(\boldsymbol{x}, t)$] として返す(結果としてベクトル計算を行わせている)。 ただし、$\boldsymbol{x} + dt \cdot\boldsymbol{fx}$ を改めて $\boldsymbol{x}$ とするために 9,10行目のように、for文を使ってリスト要素毎の代入計算 x[i] += dt * fx[i] が必要である。 このように計算できることを確かめるには(コメントアウトしている)11行目で点 $\boldsymbol{x}$ の要素をprintして確かめてみればよい。

def euler(x0, T, dt, f):
    width = len(x0)
    x = x0
    t = 0
    orbit = []
    while t <= T:
        orbit.append(list(x))
        fx = map(lambda a: a(t, *x), f)
        for i in range(width):
            x[i] += dt * fx[i]
        #print x,
        t += dt
    return orbit

可変長引数

Pythonの記号であるアスタリスク (*) は通常、乗算(掛け算)として利用される(文字列についても 'apple' * 4 は 'appleappleappleapple' と掛け算の意味として使う)。 しかし、関数定義の際の引数の並びに、アスタリスク (*) を1つ付けて可変長引数(variable arguments)の意味を持たせるように使うことができる。 アスタリスク (*) を2つ付けてキーワード可変長引数とすることもできる。 関数定義において、引数にアスタリスク (*) を1つ付ける任意の個数の引数が、1つのタプル型オブジェクト(括弧 ( と ) に挟まれたカンマ区切りの並び)の要素に格納される。

次の例を見てみよう。関数 func は形式的には x, y および arg の3つの引数を取るように見えるが、'a', 'b', 'c', 'd', 'e', 'f' の6津の引数を与えると、x に 'a' が、y に 'b' が、そして arg には残りの 'c', 'd', 'e', 'f' がセットされたことがprintで確認できる。 つまり、引数 x, y 以降に渡した任長さの引数並びが args にタプル型オブジェクトの要素として格納されたことが分かる。

>> def func(x, y, *args):
>>     print a, b, args
>> print func('a', 'b', 'c', 'd', 'e', 'f')
('a', 'b', ('c', 'd', 'e', 'f'))
>> print func('a', 'b', 'c', 'd')
('a', 'b', ('c', 'd'))
>> print func('a', 'b', 'c')
('a', 'b', ('c',))
>> print func('a', 'b')
('a', 'b', ())

Euler法でLorenz系を解いてみる

Lorenz方程式(式 (25), (26), (27))をEuler法で解いてみよう。 次のスクリプト diffeq_lorenz.py は、ベクトル場 func = [f1, f2, f3] で定まるLorenzの微分方程式を初期点 x0 = [0.1, 0.1, 0.1]、時間区間 1 、刻み幅 0.01 で(Euler法を逐次的に適用して)100個の要素からなる軌道リスト orbit を求めている。
#!/usr/bin/env python
# -*- coding: utf-8 -*-

# vector fields f1, f2, f3 of Lorenz eq as functions of t, x, y, z.
p, r, b = 10.0, 28.0, 8.0/3.0
def f1(t, x, y, z):
	return -p * x + p * y
def f2(t, x, y, z):
	return -x * z + r * x - y
def f3(t, x, y, z):
	return x * y - b * z

def euler(x0, T, dt, f):
    width = len(x0)
    x = x0
    t = 0
    orbit = []
    while t <= T:
        orbit.append(list(x))
        fx = map(lambda a: a(t, *x), f)
        #x0 = fx
        for i in range(width):
            x[i] += dt * fx[i]
            #print x,
        t += dt
    return orbit

func = [f1, f2, f3]
x0 = [0.1, 0.1, 0.1]
dt = 0.01
T = 1

#orbit = euler(x0, T, dt, func)
orbit = runge_kutta(x0, T, dt, func)
print orbit
演習: diffeq_lorenz.py を実行してみなさい。 T= 50, dt = 0.01 としたとき、 orbit の各要素(3つの要素($x$, $y$, $z$成分からなる)に関して、どれか2つだけを取り出して(3組 $(x,y), (y,z), (x,z)$ のいずれか)csv形式のファイル lorenz_euler.csv に書き出して、各点を「つない」で折れ線でプロットしてみなさい(参考:ファイルの読み書き)。

Runge-Kutta法

演習:調和振子の微分方程式で分かるように、Euler法は誤差が単純に蓄積し、真の解から急速に遠ざかってしまい、微分方程式の数値解としては実用的には使えない。

次の関数 runge_kutta(x0, T, dt, f) は、時間刻み dt の中間点を使って大幅に誤差の発生を抑える工夫がなされたRunge-Kutta法を使った関数である。 時刻 $t_0$ における$m$次元(2行目の width に次元の値が格納される)の初期点リスト x0 =[$x_{01}$,$\dots$,$x_{0m}$]から出発して、時刻 $t_0 + T$ までを指定した時間刻み dt で $N=T/dt$ 個からなる点のリスト [$\boldsymbol{x}_0$, $\boldsymbol{x}_1$, $\dots$, $\boldsymbol{x}_N$] を返す(3行目で $\boldsymbol{x}_0$ を $\boldsymbol{x}$ とおいて、while文で更新した点を orbit.append(list(x)) によって数値軌道リスト orbit に次々と追加している)。

Euler法と異なり、while文内で、 $t_n$ での値 $\boldsymbol{x}_n$ から次の時間 $t_{n+1}=t_n + dt$ での値 $\boldsymbol{x}_{n+1}$ を求める(22,23行目で与えている)ために、中間の刻み dt/2 を使った増分を計算した f1, f2, f3, f4 を使っている(計算量は4倍になる)。

def runge_kutta(x0, T, dt, f):
    width = len(x0)
    x = x0
    t = 0
    orbit = []
    while t <= T:
        orbit.append(list(x))
        x1 = x
        f1 = map(lambda a: a(t, *x1), f)
        x2 = x
        for i in range(width):
            x2[i] += dt/2 * f1[i]
        f2 = map(lambda a: a(t + dt/2, *x2), f)
        x3 = x
        for i in range(width):
            x3[i] += dt/2 * f2[i]
        f3 = map(lambda a: a(t + dt/2, *x3), f)
        x4 = x
        for i in range(width):
            x4[i] += dt/2 * f3[i]
        f4 = map(lambda a: a(t + dt/2, *x4), f)
        for i in range(width):
            x[i] += dt/6 * (f1[i] + 2 * f2[i] + 2 * f3[i] + f4[i])
        t += dt
    return orbit
演習: 先の diffeq_lorenz.py に関数 runge_kutta(x0, T, dt, f) を追加し、Runge-Kutta法で Lorenz方程式を数値計算してみなさい。 T= 50, dt = 0.01 としたとき、 orbit の各要素(3つの要素($x$, $y$, $z$成分からなる)に関して、どれか2つだけを取り出して(3組 $(x,y), (y,z), (x,z)$ のいずれか)csv形式のファイル lorenz_runge.csv に書き出して、各点を「つない」で折れ線でプロットし、Euler法による結果とを比較してみなさい(参考:ファイルの読み書き)。

モジュールGnuplot.pyを使う

利用しているコンピュータに gnuplotが使えるときには、PythonからGnuplotを使うGnuplot.pyをインストールしよう。

次は、Gnuplotモジュールをインポートしたスクリプト diffeq_lorenz.py の後半部分である。Runge-Kutta法によって軌道 orbit = runge_kutta(x0, T, dt, func) を取得して他の値に、Gnuplotでその様子をプロットするために、3つの座標成分を持つ点列データ orbit = [$\boldsymbol{x}_0$, $\boldsymbol{x}_1$, $\dots$, $\boldsymbol{x}_N$] を xorbit = [$x_0,x_1,\dots,x_N$], yorbit = [$y_0,y_1,\dots,y_N$], zorbit = [$z_0,z_1,\dots,z_N$] に成分ごとのデータに分けて(17から23行目)、これを27行目で dat = Gnuplot.Data(xorbit,yorbit,zorbit) と渡して 28行目以降で Gnuplotに描かせている。

左図は T = 50, dt =0.05 としてGnuplotでプロットした様子である。
import Gnuplot

<diffeq_lorenz.py の関数定義あれこれ>

func = [f1, f2, f3]
x0 = [0.1, 0.1, 0.1]
dt = 0.005
T = 50

####  getting numerical orbits  ####
#orbit = euler(x0, T, dt, func)
orbit = runge_kutta(x0, T, dt, func)
#print orbit

####  arranging orbit data avairable for Gnuplot  #####
xorbit = []
yorbit = []
zorbit = []
for i in range(len(orbit)):
    xorbit.append(orbit[i][0])
    yorbit.append(orbit[i][1])
    zorbit.append(orbit[i][2])
#print xorbit
#print yorbit

dat = Gnuplot.Data(xorbit, yorbit, zorbit)
gp = Gnuplot.Gnuplot()
gp('set terminal x11')
gp('set style data lines')
gp.set(xlabel="x",
 ylabel="y",zlabel="z")
gp.splot(dat)

Euler法とRunge-Kutta法の直接比較

Lorenz方程式では、原点付近から出発した解軌道は無限遠には離れないことがわかっているために、刻み時間幅 $\Delta t$ をある程度小さくして原点付近から出発させた数値解がいきなり無限編に発散してしまうようなことはない。

では、Runge-Kutta法はEuler法に比べてどれほど改良されているのかは直ちに明らかではない。 比較検討のためには、厳密解が分かっている場合について 調べてみるとよい。

演習(調和振子):

質量 $m$ の質点がバネ定数 $K=k^2 (k>0)$ のバネでつながれて直線上を運動する調和振子の微分方程式は、位置を $x$ および運動量を $p$ として次のようになる。

\begin{align*} \frac{dx}{dt} &= \frac{p}{m}\\ \frac{dp}{dt} &= -k^2 x \end{align*}

次の関数 $E(t)$ \begin{equation} E(t)=\frac{1}{2m}p(t)^2+\frac{k^2}{2}x(t)^2 \end{equation} を微分すると直ぐに確かめられるように$dE/dt=0$ つまり $E$ は一定で(エネルギー保存則)で、次式。 \begin{align*} x(t) &= \frac{\sqrt{2E}}{k} \sin\left( \frac{k}{\sqrt{m}}t + \theta_0 \right)\\ p(t) &=\sqrt{2mE}\cos\left( \frac{k}{\sqrt{m}}t + \theta_0 \right) \end{align*} は解になっていて、$(x,p)$ は楕円上を周期 $\displaystyle\frac{2\sqrt{m}}{k}\pi$ で周回する。

簡単のために、$m=1$, $k=\sqrt{2}$, $E=1$, $\theta_0=0$、つまり $x(0)=0$, $p(0)=\sqrt{2}$ を初期条件とする運動を考える(周期 $\sqrt{2}\pi$ )。

  1. Euler法とRunge-Kutta法とで求めた数値解 $(x_n,p_n)$ を(csvファイルに書き出して)それぞれプロットして、比較してみなさい。
  2. Euler法とRunge-Kutta法とで求めた数値解 $(x_n,p_n)$ の式(1)のエネルギー値 $E_n = \frac{1}{2m}p_n^2 + \frac{k^2}{2}x_n^2$ を計算し、その値を比較しなさい。

数値解の誤差

時間刻み $\Delta t$ を小さくすることが望ましいことは明かである。 正しい結果は $\Delta t\rightarrow 0$ の極限で「のみ」もたらされるのであって、微分方程式の数値解から得られ結果は真の解ではない。 しかし、$\Delta t$ を小さくすればするほど、時間経過は遅くなり計算量が増す(それでも有限の時間刻みで計算している以上、誤差の存在から逃れられない)。

また、時間刻み $\Delta t$ を一定にする積極的な理由もない。 ベクトル場 $\boldsymbol{f}(\boldsymbol{x},t)$ が険しいときには、わずかな $\Delta t$ であっても大きな変化 $\boldsymbol{x}(t+\Delta t)-\boldsymbol{x}(t)$ をもたらすため、そのような箇所では時間刻みを細かくすることが望ましいからだ。 たとえば、万有引力やCoulomb力のように距離の二乗に反比例する力が働くとき場合である。

常微分方程式の数値解法は、ここでも取り上げる中間点を考慮して精度を向上させるRunge-Kutta法など実用上広く用いられている方法があるが、誤差限界は存在している。 非線形方程式では、たとえばLorenz方程式のような単純な場合でも真の解の挙動が知られていないために、それ自体が研究対象になり得る(SmaleのMathematical Problems for the Next Century(Mathematical Intelligencer 20 (2): 7–15)の問題14)。

実際、 精度保証付き数値計算の研究は W. TuckerのA Rigorous ODE Solver and Smale’s 14th Problemによって非線形常微分方程式については一応の解決をみたが(参考:「精度保証付き数値計算の力学系への応用について」)、偏微分方程式においては依然として大きな研究テーマであり続けている。

Lorenz方程式を研究する

初期条件を与えればそれ以降の挙動は決定論的に(解の存在と一意性があれば)一意的に定まってしまう微分方程式(力学系ともいう)の研究において、 Edward N. LorenzのDeterministic Nonperiodic Flow J. Atmos. Sci., 20(11963), 130–141)は極めて大きな役割を果たした。 ほんの僅かな初期条件の違いによって、その後の挙動が極めて鋭敏に依存して全く異なってしまうという初期条件鋭敏性(sensitive dependence on initial conditions)によって、事実上挙動が予想できないというカオス(chaos)研究の先鞭をつけた。

Lorenz系をRで数値計算でも取り上げたが、数値軌道を調べてみよう。 ベクトル場のパラメータは以降 $p = 10.0$ , $r = 28.04$, $b = 8.0/3$ とする。 パラメータ依存性の研究では、$p=10, b=8/3$ を固定して、$r$ を変化させるのが通例である。

Poincareの切断面の方法 surface of section

この$1 < r$のときには、原点 $(0,0,0$ 以外に2つの不動点 $C_{\pm}$ \[ C_{\pm}=(\pm\sqrt{b(r-1)}, \pm\sqrt{b(r-1)}, r-1) \] が現れることが分かっている(z成分は2点とも$r - 1$で、同じ高さにある)。 真の軌道$(x(t), y(t), z(t))$ は平面 $z=r-1$ で定まる$xy$-平面 $\Sigma_{z=r-1}$ を「下から上に」あるいは「上から下へ」と横断(transverse)する。 たとえば、「上から下に」平面 $\Sigma_{z=r-1}$ を横断する$xy$-座標を求めてみよう。

数値軌道リスト リスト [$\boldsymbol{x}_0$, $\boldsymbol{x}_1$, $\dots$, $\boldsymbol{x}_N$] において、各要素の場所2(3番目)の値が今注目している$z$ 成分の列で $(z_0,z_1,\dots, z_N$ である。 各要素から平面 $\Sigma_{z=r-1}$ の高さを引いた列 $s=(z_0-(r-1),z_1-(r-1),\dots, z_N-(r-1))$ を考える。

$k$番目の点 $\boldsymbol{x}_k$ が平面 $\Sigma_{z=r-1}$ を「上から下に」横断するとは、 次のような $k+1$番目の点 $\boldsymbol{x}_{k+1}$ が続いているときである。 \[ \text{$\boldsymbol{x}_k >0$ and $\boldsymbol{x}_{k+1}<0$} \Leftrightarrow \text{$s_k > 0$ and $s_ks_{k+1}<0$} \] ベクトル場が滑らかで時間刻みが十分に小さければ、そのような点 $\boldsymbol{x}_k$ は平面 $\Sigma_{z=r-1}$ にごく近いと見なすことができるはずだ。 $\boldsymbol{x}_k$ の $xy$成分が、軌道と平面 $\Sigma_{z=r-1}$ とが横断する交点 $P$ であると考える。

軌道と平面 $\Sigma_{z=r-1}$ とが点 $\boldsymbol{x}_i$ で横断したとする。 平面 $\Sigma_{z=r-1}$ を $\boldsymbol{x}_i$ から出発した軌道は再び平面 $\Sigma_{z=r-1}$ を「上から下に」点 $\boldsymbol{x}_{i+1}$ で横断する。 一般に、周期軌道に横断的な平面 $\Sigma$ を考える。 平面 $\Sigma$ 上の周期軌道の近傍から出発した軌道は再び $\Sigma$ を同じ方向に横切ること、つまり平面 $\Sigma$ から $\Sigma$ への写像 \[ P: \Sigma \rightarrow \Sigma \] を定めている。 周期軌道に近傍で定義される平面 $\Sigma$ 上の写像をPoincare写像という。 周期軌道が $\Sigma$ と交差する点は写像 $P$ の不動点になっていることに注意する。

$n$次元の微分方程式の研究は、周期軌道に横断的に横切る $n-1$ 次元横断面 $\Sigma$ 上のPoincare写像の研究に帰着される。

演習: スクリプト diffeq_lorenz.py において、Runge-Kutta法を使って、T = 100、 dt = 0.01 として、平面 $\Sigma_{z=r-1}$ を上から下に横断する点を近似的に求め、その$xy$座標を csvファイル lorenz_section.csv に書き出して、プロットしてみなさい。

Lorenzプロット

E. N. LorenzはDeterministic Nonperiodic Flow J. Atmos. Sci., 20(11963), 130–141)において、今日ではLorenzプロットとして知られる力学系の対象次元をさらに減らす別の方法を提案した。

原点近くを出発した軌道は、今のパラメータの設定値では、双曲不安定点である原点を離れて暫くすると(移行期を過ぎると)、二つの$z=r-1$の高さにある2つの不動点 $C_{\pm}$ の回りを不規則的に交互に周回する(何回かを不動点 $C_+$ または $C_-$ の回りを何回転かしてから他方の不動点の回りを回り、これを交互に不規則に繰り返す)。

数値軌道リスト [$\boldsymbol{x}_0$, $\boldsymbol{x}_1$, $\dots$, $\boldsymbol{x}_N$] において、$z$ 成分の列 $(z_0,z_1,\dots, z_N)$ の隣り合う値の差 $zdif_i = z_{i+1}-z_i$に注目する。 移行期までのステップ数 transit = 500 として、transit番目以降の $z_i$ について、つぎのようにて局所極大 $\text{$z_{i-1}< z_i$ and $z_i > z_{i+1}$}$ となる$z$値の列 zextremum を求める。

import Gnuplot

....
....
....

zextremum = []
transit = 500
for i in range(transit, len(zorbit)-1):
    if zorbit[i-1] < zorbit[i] and zorbit[i] > zorbit[i+1]:
        zextremum.append(zorbit[i])
#print len(zextremum)
zpair = []
for i in range(len(zextremum)-1):
    zpair.append([zextremum[i],zextremum[i+1]])
#print 'zpair = ', len(zpair)


#### for Gnuplot
lp = Gnuplot.Data(zpair)
glp = Gnuplot.Gnuplot()
glp('set terminal x11')
glp.set(xlabel="z[i]", ylabel="z[i+1]")
glp.plot(lp)

上のスクリプトで、 zextremum =[$z_{max,1},\dots,z_{max,i-1},z_{max,i}, z_{max,i+1},\dots, z_{max, N_k}$] より、 zpair = [$\dots, (z_{max,i-1},z_{max,i}),(z_{max,i},z_{max,i+1}),\dots$] のようにして隣り合う2点を組にして $\{(z_{max,i},z_{max,i+1})\}$ をプロットして得られたのが左図である。

この図が有名なLorenzプロットである。 左図は1次元区間 $I$ から $i$自身への写像 $L$

\[ L:I \rightarrow I \] が存在することを示唆する。 もし、この写像が存在するならば、 $z_0\in I$ からテント型関数によって、$z_{n+1}=Lz_n=L^n z_0$ のように$L$を反復適用して得られる点列 $z_0,z_1,z_2,\dots$ として、Lorenz方程式の$z$方向の局所極大値が得られることになる。 この写像の存在が明らかになれば、Lorenz方程式の研究は1次元区間 $I$ 上の写像 $L$ の研究に帰着されることになる。

演習: スクリプト diffeq_lorenz.py において、上のような2点を組とする点列を csv ファイル lorenz_lpot.csv に書きだして、Lorenzプロットを描いてみなさい。