KdV方程式の解を調べる（Maximaを使ったグラフの描画）

物理の具体的な計算にMaximaを使って、物理もMaximaも同時に学んでしまいましょう。

今回はMaximaを使ってKdV方程式の解を計算してみたいと思います。

とくに、解が表す波の配位や時間発展をグラフ化・アニメ化して視覚的にわかりやすくしてみました。

グラフの作成にはMaximaに実装されている plot2d、plot3d を、

GIFアニメの作成にはgnuplotを用いてみました。

（一番最後のGIFは、Maximaで作成した複数のEPSファイルをImageMagickでGIF化しました。

というのは、gnuplotにはJacobiの楕円関数が組み込まれていないみたいだったためです）

KdV方程式
１ソリトン解
２ソリトン解
３ソリトン解
Nソリトン解
cnoidal 解

KdV方程式

KdV方程式（Korteweg-de Vries方程式）は非線形な波動を記述する方程式の一つです。

線形な波動の場合には、２つの解があると、その線形結合（重ね合わせ）もまた解になります。

しかし、非線形な波動はそのような重ね合わせの原理には従いません。

非線形波動の面白い点は、その解に粒子的な性質を持つ孤立波（ソリトン解）が出現することです。

ソリトンについての研究の歴史は長く、1834年にスコットランドの John Scott Russell が、

運河の水面を伝わる孤立波を観察したことに始まるとされています。

それから60年ほど経った1895年、オランダの D. J. Korteweg と G. de Vries により、

Russell が発見した孤立波は、こんにちKdV方程式として知られる非線形方程式

$\displaystyle{\frac{\partial u}{\partial t}+6u\frac{\partial u}{\partial x}+\frac{\partial^3 u}{\partial x^3}=0}$

によって記述されることがわかりました。

ここで $u$ は波の高さ、 $x$ は波の伝播を記述する位置座標、 $t$ は時刻に相当する変数です。

本記事では、戸田盛和『ソリトンと物理学』（SGCライブラリ49）を参考にして、

KdV方程式の解をいくつか調べてみたいと思います。

<br />

１ソリトン解

以下、略記法 $\partial_x=\partial/\partial x, ~\partial^2_x=\partial^2/\partial x^2$ などを用います。

天下り式ですが、KdV方程式の１ソリトン解（１つの孤立波からなる解）は次で与えられます：

$\displaystyle{u(x,t)=2\partial_x^2\log(1+e^{2k(x-ct)})\quad (c=4k^2)}$

これが解になっていることはKdV方程式に直接代入して

u:2*diff(log(1+exp(2*k*(x-4*k^2*t))), x, 2)$
K:diff(u,t,1)+6*u*diff(u,x,1)+diff(u,x,3)$
is(0=factor(xthru(K)));

によって確かめることができます。

また、解は $x-ct$ の関数なので、進行波を表していることがわかります。

このとき $u$ の逆数をとって整理すると

xthru(expand(ratsimp(1/u)));

より

$\displaystyle{\frac1u=\frac{e^{2k(x-ct)}+e^{-2k(x-ct)}+2}{8k^2}=\frac{\cosh\left( 2k(x-ct)\right)+1}{4k^2}}$

さらに $\cosh(2z)+1=2\cosh^2(z)$ ：

trigsimp(trigexpand(cosh(2*z)+1));

より

$\displaystyle{u=\frac{2k^2}{\cosh^2\left( k(x-ct)\right)}}$

となります。

原点 $x=0$ での波の高さの時間的推移をグラフに描画してみると（簡単のため $k=1$ としました）

plot2d (2/cosh(-4*t)^2, [t, -2, 2],[y,-0.5,2.5])$

より

（左図）からわかるように、ソリトン解はピークを持った局在した波を表しています。

右側の図は３Dプロット

plot3d (2/cosh(x-4*t)^2, [x, -1, 1],[t, -1, 1])$

によって描画したもので、孤立した波が伝わっていく様子がわかります。

また、gnuplotを起動して

set xr[-15:15]
set yr[-0.5:2.5]
set term gif animate optimize delay 2 size 480,360
set output 'movie1soliton.gif'
do for [i = -150 : 150] {
   t=i*0.025
   title(t) = sprintf("t = %.3f",t) 
   plot 2/(cosh(x-4*t)**2)
   set title title(t) font 'Times,16' 
   }
set out
set terminal wxt enhanced

と入力し、時間発展のGIFアニメを作成すると
f:id:pianofisica:20190109144437g:plain

より、孤立波がその波形を保ったまま伝わっていく様子がわかります。

ここで縦軸は波高、横軸は $x$ 座標を表しています。 $x$ 軸の正の方向に波形が伝播しています。

$\displaystyle{M_2(x,t)=\left(\begin{array}{cc} 1+\frac{1}{2k_1}e^{2k_1(x-c_1t)} & \frac{1}{k_1+k_2}e^{k_1(x-c_1t)+k_2(x-c_2t)} \\ \frac{1}{k_2+k_1}e^{k_1(x-c_1t)+k_2(x-c_2t)} & 1+\frac{1}{2k_2}e^{2k_2(x-c_2t)}\end{array} \right)}$

を導入します。ここで $c_i=4k_i^2$ です。すると２ソリトン解が

$\displaystyle{u_2(x,t)=2\partial_x^2\log\det M_2(x,t)}$

によって与えられます。実際、KdV方程式に直接代入して

M:matrix(
[1+1/(2*k1)*exp(2*k1*(x-4*k1^2*t)),1/(k1+k2)*exp(k1*(x-4*k1^2*t)+k2*(x-4*k2^2*t))],
[1/(k2+k1)*exp(k1*(x-4*k1^2*t)+k2*(x-4*k2^2*t)),1+1/(2*k2)*exp(2*k2*(x-4*k2^2*t))])$
u:2*diff(log(determinant(M)),x,2)$
K:diff(u,t,1)+6*u*diff(u,x,1)+diff(u,x,3)$
is(0=factor(factor(K)));

から、確かに解になっていることがわかります。

ここで $k_1=3/5,~k_2=1$ としたときの座標 $x=5$ における波高の時間的推移を見ると

u0:subst(1,k2,subst(3/5,k1,subst(5, x, u)))$
plot2d (u0, [t, -10, 10])$

より左図、

あるいは３Dプロット

plot3d (subst(1,k2,subst(3/5,k1, u)), [x, -10, 10],[t, -10, 10])$

（右図）から、２つの孤立波が存在し、一方が他方を追い抜いていく様子がわかります。

また、gnuplotで時間発展をGIFにしたのが以下のものです：
f:id:pianofisica:20190109121923g:plain

注目すべきは、２つの波が衝突したとき、相互作用（非線形項）のために波は互いに形を崩しますが、

そのあとは再び元の波形に戻って伝播していくという点です。

このような "粒子的な" 性質がソリトンの特徴です。

３ソリトン解

３ソリトン解は、行列

$\displaystyle{M_3(x,t)=\left(\begin{array}{ccc} 1+\frac{1}{2k_1}e^{2k_1(x-c_1t)} &\frac{1}{k_1+k_2}e^{k_1(x-c_1t)+k_2(x-c_2t)} & \frac{1}{k_1+k_3}e^{k_1(x-c_1t)+k_3(x-c_3t)} \\ \frac{1}{k_2+k_1}e^{k_2(x-c_2t)+k_1(x-c_1t)} & 1+\frac{1}{2k_2}e^{2k_2(x-c_2t)} & \frac{1}{k_2+k_3}e^{k_2(x-c_2t)+k_3(x-c_3t)} \\ \frac{1}{k_3+k_1}e^{k_3(x-c_3t)+k_1(x-c_1t)} & \frac{1}{k_3+k_2}e^{k_3(x-c_3t)+k_2(x-c_2t)} & 1+\frac{1}{2k_3}e^{2k_3(x-c_3t)}\end{array} \right)}$

（ただし $c_i=4k_i^2$ ）から

$\displaystyle{u_3(x,t)=2\partial_x^2\log\det M_3(x,t)}$

によって与えられます。実際、

M:matrix(
[1+1/(2*k1)*exp(2*k1*(x-4*k1^2*t)),
1/(k1+k2)*exp(k1*(x-4*k1^2*t)+k2*(x-4*k2^2*t)),
1/(k1+k3)*exp(k1*(x-4*k1^2*t)+k3*(x-4*k3^2*t))],
[1/(k2+k1)*exp(k2*(x-4*k2^2*t)+k1*(x-4*k1^2*t)),
1+1/(2*k2)*exp(2*k2*(x-4*k2^2*t)),
1/(k2+k3)*exp(k2*(x-4*k2^2*t)+k3*(x-4*k3^2*t))],
[1/(k3+k1)*exp(k3*(x-4*k3^2*t)+k1*(x-4*k1^2*t)),
1/(k3+k2)*exp(k3*(x-4*k3^2*t)+k2*(x-4*k2^2*t)),
1+1/(2*k3)*exp(2*k3*(x-4*k3^2*t))])$
u:2*diff(log(determinant(M)),x,2)$
K:diff(u,t,1)+6*u*diff(u,x,1)+diff(u,x,3)$
is(0=factor(factor(K)));

から、（少々計算に時間がかかりますが）確かに解になっています。

ちなみに、Maximaで

showtime:true$

とすると、計算に要した時間を計測できます。

上の計算には25秒ほどかかりました（もちろん、使用している計算機の処理能力に依存します）。

最後に、３ソリトン解を３Dプロットすると

plot3d (subst(2/5,k3, subst(4/5,k2,subst(3/5,k1, u))), [x, -20, 20],[t, -20, 20]);

より

f:id:pianofisica:20190109024054p:plain — 3-soliton 3D

と描画されて、確かに３つの孤立波が存在していることがわかります。

Nソリトン解

さて、以上の考察をまとめると、 ${N\times N}$ 行列

$\displaystyle{\left(M_N(x,t)\right)_{ij}=\delta_{ij}+\frac{\alpha_i\alpha_j}{k_i+k_j}e^{k_i(x-c_it)+k_j(x-c_jt)} \quad (c_i=4k_i^2)}$

（ここで $\alpha_i$ は任意定数です。２、３ソリトン解では簡単のため $\alpha_i=1$ としていました）

から、一般の $N$ ソリトン解が

$\displaystyle{u_N(x,t)=2\partial_x^2\log\det M_N(x,t)}$

によって構成できることがわかると思います。

cnoidal 解

孤立したソリトン解を等間隔に並べたような周期解が cnoidal 解です。

cnoidal 解は Jacobi の楕円関数を使って以下のように与えられます：

$\displaystyle{u(x,t)=2\kappa^2k^2{\rm cn}^2(k(x-ct), \kappa^2)\quad (c=-4(1-2\kappa^2)k^2)}$

ただし $\kappa$ は楕円関数の母数です。

楕円関数については以下の記事でまとめています：

pianofisica.hatenablog.com

<br />

楕円関数の性質 ${\rm cn}(z,1)=1/\cosh(z)$ から $\kappa\to 1$ で cnoidal 解は１ソリトン解に一致します。

cnoidal 解がKdV方程式の解であることは

c:-4*(1-2*kappa^2)*k^2$
u:2*kappa^2*k^2*jacobi_cn(k*(x-c*t), kappa^2)^2$
K:diff(u,t,1)+6*u*diff(u,x,1)+diff(u,x,3)$
subst(sqrt(1-kappa^2*jacobi_sn(k*(x-c*t), kappa^2)^2),jacobi_dn(k*(x-c*t), kappa^2),
subst(sqrt(1-jacobi_sn(k*(x-c*t), kappa^2)^2),jacobi_cn(k*(x-c*t), kappa^2), K))$
is(0=factor(%));

として直接KdV方程式に代入することによりわかります。ここで Jacobi の楕円関数の定義式

$\displaystyle{{\rm cn}(u,\kappa^2)=\sqrt{1-{\rm sn}^2(u,\kappa^2)} \\ {\rm dn}(u,\kappa^2)=\sqrt{1-\kappa^2\,{\rm sn}^2(u,\kappa^2)}}$

を代入しました。 $k=1$ と固定して $\kappa=~1,~1/2,~1/10$ の場合でMaximaに

plot3d (subst(1,kappa,subst(1,k, u)), [x, -1.5, 1.5],[t, -1.5, 1.5])$
plot3d (subst(1/2,kappa,subst(1,k, u)), [x, -1.5, 1.5],[t, -1.5, 1.5])$
plot3d (subst(1/10,kappa,subst(1,k, u)), [x, -1.5, 1.5],[t, -1.5, 1.5])$