KdV方程式のRによる計算 - 吾未知足唯修身爾

KdV方程式は1895年にKortewegとde Vriesによって提出された方程式であり、浅水波の減少を記述したものである。また、解の挙動が非常に面白いことでも知られる。KdV方程式は非線形偏微分方程式であるものの解が明示的に解けることでも知られるが、この記事では、数値計算法によりKdV方程式の解を求めたのち可視化し、その挙動を観察する。

なお、ここでの数値計算はRで行ったが、これは私がRを最も使い慣れていたためで（といっても触った程度だが）、特にその必然性は無い。

KdV方程式とは
KdV方程式の差分方程式
計算結果
計算結果の考察
計算コード（Rコード）

KdV方程式とは

係数の書かれ方としては色々あるが、ここではZabuskyら（1965年）が数値計算した際の式を使う。

$\dfrac{\partial u}{\partial t} + u\dfrac{\partial u}{\partial x} +\delta^{2} \dfrac{\partial^{3} u}{\partial x^{3}} = 0$

この第2項が非線形項であり、第3項は分散項と呼ばれ、波を減衰させる効果のある項である。なお、第3項が ${\partial^{2} u}/{\partial x^{2}}$ の形をした式はBurgers方程式として知られ、KdV方程式とはまた違った解の挙動を示すが、この記事では立ち入らない*1。

KdV方程式の差分方程式

続いて、数値計算用に、KdV方程式を差分方程式に書き換える。ここでもZabuskyらが使った方法（中心差分による方法、"Leap-Frog差分法"と呼ばれる）を用いる。

$u_{m}^{n+1} \\=u_{m}^{n-1} \\ -\dfrac{\Delta t}{3\Delta x}(u_{m+1}^{n}+u_{m}^{n}+u_{m-1}^{n})(u_{m+1}^{n}-u_{m-1}^{n}) \\ + \delta^{2} \dfrac{\Delta t}{3\Delta x^{3}}(u_{m+2}^{n}-2u_{m+1}^{n}+2u_{m-1}^{n}-u_{m-2}^{n})$

この式は、中心差分

$\dfrac{\partial u}{\partial t}=\dfrac{u_{m}^{n+1}-u_{m}^{n-1}}{2\Delta t}, \quad \dfrac{\partial u}{\partial x}=\dfrac{u_{m+1}^{n}-u_{m-1}^{n}}{2\Delta x}$

や、2次精度の中心差分

$\dfrac{\partial^{2} u}{\partial x^{2}}=\dfrac{u_{m+1}^{n}-2u_{m}^{n}+u_{m-1}^{n}}{\Delta x^{2}}$

の中心差分、すなわち、

$\dfrac{\partial^{3} u}{\partial x^{3}} \\=\dfrac{ \left(\dfrac{u_{m+2}^{n}-2u_{m+1}^{n}+u_{m}^{n}}{\Delta x^{2}}\right) -\left(\dfrac{u_{m}^{n}-2u_{m-1}^{n}+u_{m-2}^{n}}{\Delta x^{2}} \right)}{2\Delta x} \\= \dfrac{u_{m+2}^{n}-2u_{m+1}^{n}+2u_{m-1}^{n}-u_{m-2}^{n}}{\Delta 2x^{3}}$

から導出されるが、あわせて

$u=\dfrac{u_{m+1}^{n}+u_{m}^{n}+u_{m-1}^{n}}{3}$

なる工夫がされている。ただ、最後の点については、このような隣接する3点の平均ではなく、単に $u=u_{m}^{n}$ としても計算結果はほぼ同じである。

計算結果

$\delta=0.022$ とし、初期条件 $(1.0+\cos \pi x)$ *2および周期境界条件 $u(0,t)=u(2,t)$ のもと上記の差分方程式に従って計算する。ここで、 $dx=0.01$ 、 $dt=5.0\times10^{-4}$ とした。

この場合の解の挙動は以下の通り。最初はgifファイルを添付しようと思ったがサイズが大きかったのと、一時停止や巻戻しが出来るようyoutube経由で。

www.youtube.com

計算結果の考察

コサインカーブから始り、時間の経過と共に波が突き立ち、突き立った波がいくつかの孤立波を発生させる。

そして、大きい孤立波が小さい孤立波をすり抜けていく現象が随所に確認できる。（この孤立波はソリトンとして知られている。）

さらに時間が進むと、孤立波の位相が合ってくると数が減り、位相がずれてくるとまた増えたりするが、最後には（概ね）もとのコサインカーブに戻る。また、それを超えると再び波が突き立ち始め、いくつかの孤立波が発生し始める。

実に興味深い挙動である。

計算コード（Rコード）

こちらは計算に使ったコード。特に洗練されたものではないので参考までに。

なお、計算結果はResultという行列（x軸を行、t軸を列）に格納している。また、時間0の初期値を入れるために最初の2列（tt=1,2）を使っている。

delta=0.022
dx=0.01
dt=5*10^-4

#計算範囲の指定(x=0~2,t=0~10)
xx=2/dx
tt=10/dt
result=matrix(0, nrow=xx, ncol=tt+2)

#初期条件の指定(1+cosπx)
for (m in 1:xx){
        result[m,1]=1+cos(pi*dx*(m-1))
        result[m,2]=1+cos(pi*dx*(m-1))
}

#周期境界条件(U(0,t)=U(2,t))の下で計算
for (n in 2:(tt+1)){
        result[1,n+1]=result[1,n-1]-dt/(3*dx)*(result[2,n]+result[1,n]+result[xx,n])*(result[2,n]-result[xx,n])-delta^2*dt/dx^3*(result[3,n]-2*result[2,n]+2*result[xx,n]-result[xx-1,n])
        result[2,n+1]=result[2,n-1]-dt/(3*dx)*(result[3,n]+result[2,n]+result[1,n])*(result[3,n]-result[1,n])-delta^2*dt/dx^3*(result[4,n]-2*result[3,n]+2*result[1,n]-result[xx,n])
        for (m in 3:(xx-2)){
                result[m,n+1]=result[m,n-1]-dt/(3*dx)*(result[m+1,n]+result[m,n]+result[m-1,n])*(result[m+1,n]-result[m-1,n])-delta^2*dt/dx^3*(result[m+2,n]-2*result[m+1,n]+2*result[m-1,n]-result[m-2,n])
        }
        result[xx-1,n+1]=result[xx-1,n-1]-dt/(3*dx)*(result[xx,n]+result[xx-1,n]+result[xx-2,n])*(result[xx,n]-result[xx-2,n])-delta^2*dt/dx^3*(result[1,n]-2*result[xx,n]+2*result[xx-2,n]-result[xx-3,n])
        result[xx,n+1]=result[xx,n-1]-dt/(3*dx)*(result[1,n]+result[xx,n]+result[xx-1,n])*(result[1,n]-result[xx-1,n])-delta^2*dt/dx^3*(result[2,n]-2*result[1,n]+2*result[xx-1,n]-result[xx-2,n])
}

続いて計算結果を描画する。発想としてはggplotで時間毎に面グラフを書き、それを連続表示させるというもの。スマートなコードではないため結構無駄も多いと思われる。ともかく、まずはggplot用にデータ系列を整えつつdataframe化する。

#描画用データの整形
x_index=rep(0, times=xx)
for (m in 1:xx){
        x_index[m]=dx*(m-1)
}
x_index2=c(x_index)
t_index2=rep(0, times=xx)
result2=c(result[,2])
#全データを描画すると時間かかりすぎるので間引く
for(k in 1:1000){
        x_index2=c(x_index2, x_index)
        t_index2=c(t_index2, rep((tt/1000*k)*dt, times=xx))
        result2=c(result2, result[,2+tt/1000*k])
}
dat2=data.frame(result2, x_index2, t_index2)
dat3=dat2[dat2$t_index2!=0,] #下記グラフではt=0値が使えないため

次にggplotで時間毎に面グラフを書き（波なので面グラフを採用）、アニメーションを保存する。

#描画（計算結果のアニメーション表示）
library(ggplot2)
library(gganimate)
library(transformr)
library(magrittr)
aplot=ggplot(dat3, aes(x=dat3$x_index2, y=dat3$result2))  +
        geom_area(fill="skyblue")+
        labs(x="",y="")+
        theme_classic()+
        xlim(0, 2) + ylim(0, 4)+
        labs(title = 'time:{frame_time}', x = 'x', y = 'u(x,t)')+ 
        transition_time(dat3$t_index2)
a=animate(aplot, duration=100, fps=10, width=1280, height=720)
anim_save("kdv.wmv", a)

かつてKdV方程式をCで計算したことがありましたが、今回Rで計算してみました。間違い等ありましたらお手数ですがコメント等いただけますと幸甚です。

*1:Burgers方程式の記事はこちら

*2:計算に際しては簡便に $u(x,0)=(1.0+\cos \pi x)$ および $u(x,-\Delta t)=(1.0+\cos \pi x)$ としている。