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

KdV方程式を解く

KdV方程式の厳密解

Korteweg–de Vries (KdV)方程式は、965年のZabusky and Kruskalに数値的発見をもたらし、さらに1967年に Gardner, Greene, Kruskal, and Miuraによる僅か3ページの論文Method for Solving the Korteweg-deVries Equationにおいて、ソリトン解を有する初期値問題が逆散乱法によって解けることが示されるなど、以降の数々の数学的発見の先駆となった最も深く研究されてきた非線形偏微分方程式である。 KdV方程式はしばしば次のように表される。 \[ \frac{\partial u}{\partial t}+6u\frac{\partial u}{\partial x}+\frac{\partial^3 u}{\partial x^3}=0 \]

KdV方程式の表記の仕方はさまざまある。 $u\rightarrow -u$ の変換によって、$u_t-6uu_x+u_{xxx}=0$ や $u\rightarrow u/6$ によって $u_t+uu_x+u_{xxx}=0$ になる。 以下のZabusky and Kruskalの流儀もその1つである。

KdV方程式の特殊解

波が形を変えないで伝搬するとき、進行波（traveling wave）または定常波（stationary wave）という。 進行波は、ある関数 $f$ を使って \[ u(t,x)=f(x-ct) \] と表される。 $c$ は波の伝搬速度である。 KdV方程式の特解として、進行波解 \[ u(t,x)=\frac{c}{2}\mathrm{sech}^2 \frac{\sqrt{c}}{2}(x-ct+\delta) \] を持つことは直接代入によって確認することができる。 この解をソリトン解と呼ぶ（ただし、ソリトン解を有する他の方程式が $\mathrm{sech^2}$ の関数形を解として持つわけではない。その鋭く立ち上がって（微係数を含む）孤立した波の形状を指す意味で使う）。 ここで、$\mathrm{sech}$ は双曲線余弦関数 $\cosh$ の逆数 \[ \mathrm{sech}\, x = \frac{1}{\cosh x} =\frac{2}{e^x + e^{-x}} \] で、 \[ c=4\kappa^2 \] の関係にある。

In[1]:=EqKdV[u_, t_, x_] := 
 D[u[t, x], t] + 6 u[t, x] D[u[t, x], x] + D[u[t, x], {x, 3}] == 0
In[2]:=stationary[t, x] := c /2  Sech[Sqrt[c]/2 (x - c t + \[Delta])]^2
In[3]:= Simplify[EqKdV[stationary, t, x]]

Out[3]= True

KdV方程式の数値的研究

Zabusky and Kruskalは1965年の記念碑的論文 Interaction of "Solitons" in a Collisionless Plasma and the Recurrence of Initial States（Phys. Rev. Lett. 15, 240: ここからでも閲覧できるようだ）において、 KdVがソリトン解を有することを強く示唆する数値研究を行い、Gardner, Greene, Kruskal, and Miuraによる理論的研究を先導した。

Zabusky and KruskalはKdV方程式 \[ u_t + uu_x +\delta^2u_{xxx}=0 \] を、周期境界条件 \[ u(t,x)= u(t,x+L) \] のもとで（$L=2$）、時刻 $t=0$ での初期値問題 \[ u(0,x) =\cos \pi x \] を数値的に解析し、ソリトン（soliton: 孤立波）現象を突き止めた。 $u(t,x)$ を時刻 $t$ と位置 $x$ における波高とみなしたとき、 余弦関数を波高の初期値と与えた波が時間と共にソリトン（soliton: 孤立波）に別れて、互いに粒子のように衝突するというのである（是非、論文のFig.1 を眺めて欲しい）。 今日では、KdV方程式はN-soliton解を持つことが分かっている。

ソリトン（soliton）とは -on が付いているように、互いの衝突によっても、あたかも孤立した（solitary）粒子であるかのようにすり抜けて衝突後に元の波形を回復することが命名された特別な意味を持つ名称で、solitary waveともいう。

Zabusky-Kruskalの方法

偏微分方程式の数値解法は、誤差評価において明らかに常微分方程式以上に難しく、計算量の観点からも問題があり、今日においても主要な数値計算上の研究テーマである。

$u(t,x)=u(j\Delta t, i\Delta x)$ は \begin{align*} u\left((j+1)\Delta t,i\Delta x\right) &=\\ &=\mathcal{F} \left[\left\{ u\left((j-1)\Delta t, i\Delta x\right), u\left(j\Delta t, (i+2)\Delta x\right), u\left(j\Delta t, (i+1)\Delta x\right),\\ u\left(j\Delta t, i\Delta x\right), u\left(j\Delta t, (i-1)\Delta x\right), u\left(jd\Delta t, (i-2)\Delta x\right) \right\}\right] \end{align*} から逐次的に求められると考える。

ここでは、Zabusky and Kruskalの方法に準拠して結果だけを紹介する。 Interaction of "Solitons" in a Collisionless Plasma and the Recurrence of Initial Statesでは、$u(t,x)=u(j\Delta t, i\Delta x)= u^j_i$, $i=0,1,\dots, 2N-1$とし、周期境界条件 $u(t,x=0)=u(t,x=2)$ のもとで次のような差分スキームを使った（論文末尾の注(6)）。 \begin{align*} u^{j+1}_i &= u^{j-1}_i-\frac{\Delta t}{3\Delta x}\left(u^j_{i+1}+u^j_i+u^j_{i-1}\right)\left(u^j_{i+1}-u^j_{i-1}\right) -\delta^2 \frac{\Delta t}{\Delta x^3}\left(u^j_{i+2}-2u^j_{i+1}+2u^j_{i-1}-u^j_{i-2}\right)\\ u^j_i &= u^j_{i+2N}\qquad \text{periodic boundary condition} \end{align*} Zabusky and Kruskalでは \[ \delta=0.022,\quad \Delta t=\Delta x=1/N \] を採用している（$N=128$ 程度）。

ここでは、KdV方程式を \[ u_t + \epsilon\, uu_x +\mu\, u_{xxx}=0 \] とするとき、$u(jdt,idx)$ に関する差分スキームにおける位置 $x=i \Delta x$ の添字範囲を \[ i=0,1,\dots , N \] とすることにして（Z-Kでは i=0,1,,,.2N-1 としている）、時間発展については Zabusky and Kruskalにならい、 $i=0,1,\dots, N-1$ として 次の差分スキーム (ZK)とする。 \begin{align*} u^{j+1}_i &= u^{j-1}_i-\frac{\epsilon\Delta t}{3\Delta x}\left(u^j_{i+1}+u^j_i+u^j_{i-1}\right)\left(u^j_{i+1}-u^j_{i-1}\right) -\frac{\mu\Delta t}{\Delta x^3}\left(u^j_{i+2}-2u^j_{i+1}+2u^j_{i-1}-u^j_{i-2}\right), \tag{ZK}\\ u^j_i &= u^j_{i+N}\qquad \text{periodic boundary condition} \end{align*} ただし、初期時間 $j=0$（$t=0\Delta t$） から $j=1$（$t=\Delta t$） については、次の時刻1経過後を求めるスキーム (*) を使うことにする。 \begin{equation} u^1_i = u^0_i-\frac{\epsilon\Delta t}{6\Delta x}\left(u^0_{i+1}+u^0_i+u^0_{i-1}\right)\left(u^0_{i+1}-u^0_{i-1}\right) -\frac{\mu\Delta t}{2\Delta x^3}\left(u^0_{i+2}-2u^0_{i+1}+2u^0_{i-1}-u^0_{i-2}\right) \tag{*} \end{equation}

Pythonでの実装

パラメータの設定

$\epsilon = 0.2=$e, $\mu = 0.1=$mu, N$= 256$, $\Delta x= 0.4x=$dx, $\Delta t = 0.1=$dt とし、時間間隔をたとえば、T = 300 程度として、Python での実装を考えよう。

計算データの保持の仕方を考える

差分スキームを子細に検討するとわかるように、時刻 $j+1$ を計算するために、時刻 $j$ と時刻 $j-1$ のデータを使っている。 時刻 $t=j\Delta t$ での各点 $x=i\Delta x$ の値 $u^j_i$ （$i=0,\dots, N-1$）を長さ N のリスト u[0], u[1], .... として保持し、これを差分スキームに従って、時刻 $(j+1)\Delta t$におけるの各点 $x=i\Delta x$ での値 $u^{j+1}_i$ として逐次的に求めているわけだ。 \begin{align*} u[2] &\leftarrow u[0], u[1]\\ u[3] &\leftarrow u[1], u[2]\\ \vdots\\ u[n+1] &\leftarrow u[n-1], u[n]\\ \vdots \end{align*} このとき、全ての時刻 $j=0,1,\dots$ についてのリストを保持することはメモリ上の困難が予想される。 ここでは u[0], u[1], u[2] のリストだけを使って（必要な空間計算量は $N\times 3$ である）、最新の u[0] だけをファイルに書き出すようにして計算することを考える。 次は、準備のために、$N$個からなる空リスト u[0], u[1], u[2] を準備する。

u =[[[] for x in range(N)] for y in range(3)]

初期条件と時刻1経過後の計算

次にリスト u[0] に位置 $x=i\Delta x$ における初期条件を代入する。

for i in range(N):
    u[0][i]= math.cos(2.0 * math.pi * i / N)

時刻1経過後のリスト u[1] の要素をスキーム (*) にしたがって、次のように計算する。 ただし、周期境界条件 $u^j_i=u^j_{i+N}$ を考慮してforでは i =2,..., N-3 まで動かすとし、for文の外側で、i=0,1, N-2, N-1 を別途計算していることに注意しよう。

import math

for i in range(2,N - 2):
    u[1][i] = u[0][i] - e * dt/(6.0 * dx) * (u[0][i+1] + u[0][i] + u[0][i-1]) * (u[0][i+1] - u[0][i-1]) - mu * dt/(2.0 * (dx ** 3)) * (u[0][i+2] - 2.0 * u[0][i+1] + 2.0 * u[0][i-1] - u[0][i-2])
# taking account of boundary condition
u[1][0] = u[0][0] - e * dt/(6.0 * dx) * (u[0][1] + u[0][0] + u[0][N-1]) * (u[0][1] - u[0][N-1]) - mu * dt/(2.0 * (dx ** 3)) * (u[0][2] - 2.0 * u[0][1] + 2.0 * u[0][N-1] - u[0][N-2])

u[1][1] = u[0][1] - e * dt/(6.0 * dx) * (u[0][2] + u[0][1] + u[0][0]) * (u[0][2] - u[0][0]) - mu * dt / (2.0 * (dx ** 3)) * (u[0][3]- 2.0 * u[0][2] + 2.0 * u[0][0] - u[0][N-1])

u[1][N-2] = u[0][N-2] - e * dt / (6.0 * dx) * (u[0][N-1] + u[0][N-2] + u[0][N-3]) * (u[0][N-1] - u[0][N-3]) - mu * dt/(2.0 * (dx ** 3)) * (u[0][0] - 2.0 * u[0][N-1] + 2.0 * u[0][N-3] - u[0][N-4])

u[1][N-1] = u[0][N-1] - e * dt/ (6.0 * dx) * (u[0][0] + u[0][N-1] + u[0][N-2]) * (u[0][0] - u[0][N-2]) - mu * dt/(2.0 * (dx ** 3)) * (u[0][1] - 2.0 * u[0][0] + 2.0 * u[0][N-2] - u[0][N-3])

時間発展

差分スキーム(ZK)と1経過後を与えるスキーム (*) を使って、次のようにして時間発展を計算する。 2行目で$N$ 個の要素からなるリスト u[0] 内容を「1行でファイルに書き出す」準備をし、while文内の7,8行目で計算回数が 15の倍数のときの u[0] の内容を書き出している（どのように値を書き出すかは考察する必要がある）。 23,24,25行目で、リスト u[0], u[1] を u[2] を使って更新している。

# ファイルに書き出す準備
fh = open('ファイル名', 'w')

t = 0
iteration = 0 # 計算回数（時刻 j）
while t <= T:
    if iteration % 15 == 0: # 15回に1つずつ書き出し
        # 数値リスト u[0] を一行でファイルに書き出す
        #for i in range(N):
        	#fh.write(str(i)+' ' + str(u[0][i]) + '\n' )
        #fh.write('\n')
        

    for i in range(N - 2):
        u[2][i] = u[0][i] - e * dt /(3.0 * dx) * (u[1][i+1] + u[1][i] + u[1][i-1]) * (u[1][i+1] - u[1][i-1]) - mu * dt / (dx ** 3)*(u[1][i+2] - 2.0 * u[1][i+1] + 2.0 * u[1][i-1] - u[1][i-2])
    # near boundary
    u[2][0] = u[0][0] - e * dt /(3.0 * dx) * (u[1][1] + u[1][0] + u[1][N-1]) * (u[1][1] - u[1][N-1]) -  mu * dt / (dx ** 3)*(u[1][2] - 2.0 * u[1][1] + 2.0 * u[1][N-1] - u[1][N-2])
    
    u[2][1] = u[0][1] - e * dt /(3.0 * dx) * (u[1][2] + u[1][1] + u[1][0]) * (u[1][2] - u[1][0]) -  mu * dt / (dx ** 3)*(u[1][3] - 2.0 * u[1][2] + 2.0 * u[1][0] - u[1][N-1])
    
    u[2][N-2] = u[0][N-2] - e * dt /(3.0 * dx) * (u[1][N-1] + u[1][N-2] + u[1][N-3]) * (u[1][N-1] - u[1][N-3]) -  mu * dt / (dx ** 3)*(u[1][0] - 2.0 * u[1][N-1] + 2.0 * u[1][N-3] - u[1][N-4])
    
    u[2][N-1] = u[0][N-1] - e * dt /(3.0 * dx) * (u[1][0] + u[1][N-1] + u[1][N-2]) * (u[1][0]- u[1][N-3]) -  mu * dt / (dx ** 3)*(u[1][1] - 2.0 * u[1][0] + 2.0 * u[1][N-2] - u[1][N-3])

    # we have only u[0][x] qnd u[1][x] for saving memory
    # update u[0][x] and u[1][x] for next computation u(x,2) and u(x,3)
    for i in range(N):
            u[0][i]=u[1][i]
            u[1][i]=u[2][i]
    t += dt
    iteration += 1

fh.close()

ソリトンの発生

周期境界条件 $u(t,x)=u(t,x+2)$ のもとで、$t=0$ での初期条件 $u(0,x)=\cos \pi x$ を与えたときのKdV方程式の数値結果を表示してみよう（N=256, dx=0.4, dt=0.1）。 表示の仕方については吟味すべきである。

以下の図（クリックで拡大）は、単純な表示としてある瞬間における波ショットとして、ファイルに書き出した各行の$N$個（この例では256個）の波高データの幾つかをそれぞれ取り出してを、MathematicaのListPlot（Joined->True オプションを付けて線で結んだ）でプロットしたものである。 600ステップくらいから、ソリトンに成長する立ち上がりが発生し、900ステップには2,3個のソリトンが現れる様子が分かる。 1200ステップになるとソリトン列となって進行していく様子（周期境界条件によって $x=0$ は $x=2$ とはつながっていると見なして観察する）。

今の場合、時間幅を $T$ として、$T/dt$ 回の時間刻みリストデータを 15回に1回ずつファイルに書き出しているため、ファイルは約 $T/(15dt)$ 行からなる。

ソリトンの発生をわかりやすく表示する