【例題】KdV方程式¶

警告

このドキュメントは古いバージョンです。最新バージョンはこちら。

方程式 ¶

以下の非線型偏微分方程式はKorteweg-de Vries (KdV)方程式と呼ばれます：

(1)¶ $\frac{\partial u(x,t)}{\partial t} + 6 u(x,t) \frac{\partial u(x,t)}{\partial x} + \frac{\partial^3 u(x,t)}{\partial x^3} = 0$

$u(x,t)$ は実数とします。各項の係数や符号は変数変換で変わるので、重要ではありません。慣例的に第2項の係数を6とする文献が多いので、ここではそれに従っています。この方程式は 孤立波 （ソリトン と呼ばれる）を解として持つことが知られています。この方程式を数値計算によって解いて、孤立波を確認してみましょう。

方程式の解説 ¶

KdV方程式を含む非線形方程式については文献 [4] をお勧めします。以下はその概要です（かなり大雑把ですが）。

KdV方程式は水深の浅い水面を進む波（浅水波）を記述する方程式として、1895年にKortewegとde Vriesによって定式化されました。(1) 式の第2項（非線形項）は波の上部ほど速度が速くなって波形が前のめりになる効果（波の突っ立ち）を表します。第3項は波を広げる役割があり、分散項と呼ばれます。非線形項の波形を崩す効果と分散項の波形を広げる効果がうまくつり合って、KdV方程式は非線形方程式であるにもかかわらず、波形の形を保ったまま進む進行波を解として持ちます。

KdV方程式は非線型方程式なので、解の重ね合わせが成り立ちません。しかし、2つの進行波が衝突して互いに離れた後、それぞれの波形を保ちます。このように、この進行波は粒子的な振る舞いをすることから、solitary waveまたはsolitonと呼ばれます。

KdV方程式は1950年代以降のフェルミ・パスタ・ウラムの問題 (FPU問題) に関連した研究で脚光を浴びました。フェルミらは、統計力学におけるエルゴード性の検証のために、非線形性を考慮した連成振動模型を考えました。通常の（線形の）連成振動の場合、各振動モードが独立に運動するため、そのエネルギーは保存されますが、非線型性を考慮した場合にはモード間に結合が生じ、モード間でエネルギーのやり取りが起こります。そのため、充分時間が経てば、初期状態に依らず、全てのモードに等しくエネルギーが分配された平衡状態に達する（熱化する）と予想できます。しかし、数値計算の結果、平衡状態への熱化は起きず、一定時間後に初期状態に戻る 再起現象 が観測されました。この現象は、非線形方程式における孤立波の存在を示唆します。実は、FPU問題の非線形連成振動模型は、ある極限でKdV方程式に一致し、KdV方程式と同じ孤立波を解として持ちます [4]。

数値解法 ¶

偏微分方程式の解法に従って、空間座標 $x$ を離散化し、 $u(x,t)$ を $\bm{u}(t)$ とベクトル表示します。ここで $[\bm{u}(t)]_i \equiv u(x_i,t)$ です。すると、方程式 (1) は

(2)¶ $\frac{d \bm{u}(t)}{dt} = -6 \bm{u}(t) \circ D_{x} \bm{u}(t) - D_{xxx} \bm{u}(t)$

と表されます。ここで、 $\circ$ はベクトルの成分同士の掛け算を表す記号（アダマール積）、 $D_{x}$ と $D_{xxx}$ は差分を表す行列です。差分の詳細は差分を参照してください。この方程式は変数 $t$ に関する常微分方程式なので、常微分方程式の解法を使って解くことができます。

実装 ¶

ソースコード kdv_plain.py

まずは差分を表す行列 $D_x$ と $D_{xxx}$ を作ります。簡単のため1次の中心差分公式を適用すると、3階微分までの差分行列は以下のように実装できます：

ソースコード differential.py

from scipy.sparse import csr_matrix as sparse_matrix

def make_differential_ops(nx, dx):

    # operators (matrix) that shift vector components (periodic boundary condition)
    f0 = np.identity(nx, dtype=int)  # f_{i}
    f1 = np.roll(f0, 1, axis=1)  # f_{i+1}
    f2 = np.roll(f0, 2, axis=1)  # f_{i+2}
    f_1 = f1.transpose()  # f_{i-1}
    f_2 = f2.transpose()  # f_{i-2}

    # (f_{i+1} - f_{i-1}) / (2 dx)
    deriv1 = sparse_matrix(f1 - f_1) / (2.0 * dx)

    # (f_{i+1} - 2f_{i} + f_{i-1}) / (dx^2)
    deriv2 = sparse_matrix(f1 - 2.0 * f0 + f_1) / np.square(dx)

    # (f_{i+2} - 2f_{i+1} + 2f_{i-1} - f_{i-2}) / (2 dx^3)
    deriv3 = sparse_matrix(f2 - 2.0 * f1 + 2.0 * f_1 - f_2) / (2.0 * np.power(dx, 3))

    return deriv1, deriv2, deriv3

2階微分は今回は使用しませんが、ついでに実装しています。この関数は最初に1度呼び出すだけなので、計算時間にはほとんど影響しません。

今回作る行列は規則的なものなので、単位行列（numpy.identity で生成）から基本的な操作を組み合わせることで、for文を使わずに書けます。ここではまず、 numpy.roll を使って単位行列を右方向(axis=1 で指定)に1つまたは2つずらした行列（f1 および f2）を作っています。このとき右にはみ出した要素は左に戻ります。これは周期境界条件に対応しています。

最終的に得られる行列は疎行列なので scipy.sparse モジュールを利用します。疎行列クラスとしていくつか選択肢がありますが、ここでは scipy.sparse.csr_matrix (Compressed Sparse Row format)を使用しています。右からベクトルを書ける演算が中心なので、行列を行ベクトルとして持っていた方が若干効率がよさそうという理由です（大差はないと思います）。

make_differential_ops で作った行列を使うと、空間微分を差分で置き換えたKdV方程式 (2) の右辺は以下のようにシンプルに書けます：

def f_kdv(u, t, df1, df3):
    u_x = df1.dot(u)
    u_xxx = df3.dot(u)
    return -6.0 * u * u_x - u_xxx

この関数はベクトル $\bm{u}(t)$ を受け取って、ベクトル $d\bm{u}(t)/dt$ を返します。引数 t は関数内で使用しませんが、ソルバーの仕様に合わせるために必要です。その他の引数 df1 と df3 は (2) 式右辺の計算に必要なので、パラメーターとして関数に与えます。

あとはこの関数を常微分方程式を解く関数 solve_ode に渡せば解くことができます。

# differential operators
op_df1, _, op_df3 = make_differential_ops(nx, dx)

print("Solving equation...")
algo = "runge_kutta"
args = (op_df1, op_df3)
u_tx = solve_ode(f_kdv, u0, t_array, args=args, algo=algo)

2行目左辺のアンダーバー '_' は、関数の戻り値を破棄するという意味です。2階微分は使用しないので2番目の戻り値を捨てています。5行目では、計算時間の節約のためにルンゲ・クッタ法を指定しています。odeint を使いたい場合は、algo = "odeint" とします。

ソースコードの最後の部分で、結果をアニメーションGIFとして出力しています。

t_skip = 100
u_tx_skip = u_tx[::t_skip, :]

u_tx という変数が $u_i(t)$ の全データを2次元配列として保持していますが、これをそのままアニメーションにすると動きがとても遅くなってしまいます。u_tx から100ステップごと（t_skip で指定）にデータを切り出してからアニメーションを作成します。

u_tx は100MBを超える大きなデータなので、数値データとして全て保存することはしていません。もし、データを保存したい場合には、numpy.savetxt などを使ってデータを保存してください。

結果 ¶

下図は初期条件をsin関数とした場合の結果です。

いくつものピークに分かれて、それらが波のように個別に進行していく様子が分かります。これを孤立波と呼びます。

不安定性 ¶

偏微分方程式の場合には、アルゴリズム特有の不安定性に注意が必要です。これは「誤差の爆発」と言われたりします。変数の離散化を適切に行わないと、たとえプログラムが正しくても、正しい結果が出ません。一見するとプログラムのバグのように見えるので、これを知らないと、存在しないバグを探し続けることになってしまいます。

KdV方程式は非線型方程式なので、標準問題のように解析的にきれいな形には書けませんが、主要部分のみを評価すると

$\frac{u_i(t+\Delta t)}{u_i(t)} \sim \frac{\Delta t}{\Delta x^3}$

となります。したがって、精度を上げようとして空間メッシュを増やした場合（ $\Delta x$ を小さくした場合）、それに応じて時間発展を少しづつ（ $\Delta t$ を小さく）進めなければなりません。

【発展】クラスを使った実装 ¶

さて、ソリトンの性質を調べるために、初期条件を変えて計算したいとします。先ほどのスクリプトをコピーして必要個所を書き換えれば計算可能です。しかし、この方法では、共通箇所を変更した場合に（バグ修正をしたり、コードを改良したり）、コピーしたファイルの全てを直さなければならず、とても面倒ですし必ずミスが起こります。プログラミングでは、コードの重複は極力避けなければいけません。そこで、初期条件に依らない共通の部分は1つのファイルにまとめておき、初期条件やパラメーターを別ファイルに持つようにスクリプトを改良することを考えます。

ここではクラスを使って実装してみます。クラスの中身を見る前に、クラスを使った場合にmain部分がどのように書けるかを見てみます。

ソースコード kdv_sin.py

from __future__ import print_function

import numpy
import math
from kdv import KDV


if __name__ == '__main__':

    nx = 1000
    x_min = 0
    x_max = 100.0

    nt = 10001
    t_max = 10.0

    # "runge_kutta" (default), "euler", or "odeint"
    algo = "runge_kutta"

    # initialize
    kdv = KDV(nx, x_min, x_max, nt, t_max)

    # set initial condition
    x = kdv.x
    u0 = np.sin(x * (2.0 * math.pi / x_max))
    kdv.u0 = u0

    # solve equation
    kdv.solve(algo)

    # animation
    t_skip = 100
    kdv.save_animation(t_skip, ylim=[-1.5, 3.0], filename="kdv_sin.gif")

計算の流れは以下の通りです：

kdv = KDV(nx, x_min, x_max, nt, t_max) でクラスを初期化し、そのオブジェクトを kdv として保持。
kdv.u0 = u0 で初期条件をセット。それに先立ち、x = kdv.x で空間メッシュを受け取る。
kdv.solve() で計算を実行。このとき引数で使用するアルゴリズムを指定する。
kdv.save_animation() で動画を作成。

先の例と比べて全体の見通しが良いことが分かると思います。なお、同様のことは関数を使っても実装できますが、その場合は関数同士で明示的にデータを受け渡しする必要があります。クラスの場合は、必要なデータをkdvオブジェクトが保持しているので、使う側はデータの受け渡しを意識しなくて良いという利点があります。

KDVクラスは次のファイルに定義されています。

ソースコード kdv.py

説明は省きますが、"self."で始まる変数がクラス変数で、この変数はkdvオブジェクトが存在する間は保持されます。それにより関数（クラスの場合はメソッドと呼ぶ）の間でデータのやり取りを行うことができます。