RでKdV方程式を数値解析でプロットしてみる

浅水波を表すKdV方程式は、今では可積分系として閉じた解が知られている非線形偏微分方程式ですが、これは意外に最近知られるようになった事で、二十世紀の中ごろにザブスキーとクルスカルと言う物理学者がシミュレーションをしてソリトンと言う孤立波を見つけ、そこから解析が大きく進むようになったそうです（武部(2008)）。当時を偲んでRで数値計算してみましょう。

1. KdV方程式の形状

まずは式の確認です。

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

$t$ が時間、 $x$ が位置、 $u(t,x)$ が波の高さになります。
このKdV方程式は、例えば波動方程式と比べると厄介な形状をしています。三階微分が中にあり、数値演算のお決まりのパターンである常微分方程式にならないのです。

2. 数値演算が出来る形に変形する

どうしたら良いのかな・・・と思っていたら、Numerical Solution of the KdVと言うページに全てのノウハウが書いてありました。
まずは第二項を整理します。合成関数の微分を思い出すと何をしているのか分かります。

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

ここでフーリエ変換を思い出します。Rやmatlabで使われる式とは周期と積分範囲が異なるので注意してください。

$F\bigg(u(k) \bigg ) = \hat{u(k)} = \int^{\infty}_{-\infty} u(x) e^{-ikx} dx$

フーリエ変換をかけます。公式により虚数 $i$ と $k$ で $x$ の偏微分を無くす事ができます。

$\frac{\partial}{\partial t} \hat{ u(k) } + 3ik \hat{ u(k)^2 } - ik^3 \hat{ u(k) } = 0$

この式はスプリット・ステップ法（the split step method）で計算する事ができます。

$\frac{\partial}{\partial t} \hat{ u(k) } = - 3ik \hat{ u(k)^2 } + ik^3 \hat{ u(k) }$

オイラー法で $\Delta t$ を用いて書くと以下のようになります。

$\hat{ u_1(k, t + \Delta t) } = \hat{ u(k, t) } e^{ik^3 \Delta t}$

$\partial \hat{u}/\partial t \cdot 1/\hat{u}=ik^3$ の形にして、両辺を $t$ で積分し、対数を消して $u(k, t)=C e^{ik^3 t}$ の形に整理してから、 $u(k, t + \delta t)$ を計算していますね。

$\hat{ u(k, t + \Delta t) } = \hat{ u_1(k, t + \Delta t) } - 3ik \Delta t \hat{ u_1(k, t + \Delta t)^2 }$

上式の右辺第二項のフーリエ変換周りは注意してください。

$\hat{ u(k, t + \Delta t) } = \hat{ u_1(k, t + \Delta t) } - 3ik \Delta t F [ F^{-1} [ \hat{ u_1(k, t + \Delta t)^2 }$ ] ]

フーリエ変換、逆フーリエ変換、フーリエ変換の三段構えです。

3. ソースコード

matlabのコードをRに翻訳した感じですが、意外に短いコードで計算する事ができます。計算速度は・・・遅いですね。なお精度はイマイチなので、ルンゲ＝クッタ法などに書き直した方が良いかも知れません。

# 補助的に使う関数
cosh <- function(x){ (exp(x)+exp(-x))/2 }
sech <- function(x){ 1/cosh(x) }
phi <- function(x, t, k=1){
  2*k^2*sech(k*(x-4*k^2*t))^2
}
ifft <- function(x){ fft(x, inverse=TRUE)/length(x) }

# 初期値
N <- 256
x <- seq(-10, 10, length.out=N)

# 波の初期値
# ソリトン波を指定しているが、他のモノでも良い
u1 <- phi(x, -1)

# 変化幅
delta_x <- x[2] - x[1]
delta_t <- 10/N^2 # エラーが出る場合は0.4/N^2で
# フーリエ変換の定義による周期の違いを2*piで調整
delta_k <- 2*pi/(N*delta_x)

# 離散フーリエ変換ではkの値が不連続になるので注意
# 参考: http://www.made.gifu-u.ac.jp/~tanaka/LectureNote/numerical_analysis.pdf P.109
k <- c(seq(0, N/2*delta_k, delta_k), seq(-(N/2-1)*delta_k, -delta_k, delta_k))

# t=1までnmax回波を動かす
tmax <- 1
nmax <- round(tmax/delta_t)
U = fft(u1)
for(n in 1:nmax){
  # first we solve the linear part
  U = U*exp(1i*k^3*delta_t)
  # then we solve the non linear part
  U = U - 3i*k*delta_t*(fft(Re(ifft(U))^2))
}
# 変化後の波を得る
u2 <- Re(ifft(U))

# 元と変化後の波を描画する
xlim = c(0, N)
ylim = c(0, max(u1, u2))
plot(u1, xlim=xlim, ylim=ylim, type="l", lty="dotted", xlab="x", ylab="u(t, x)", main="KdV方程式のプロット")
par(new=T)
plot(u2, xlim=xlim, ylim=ylim, type="l", xlab="", ylab="")
legend("topleft", lty=c(3, 1), col=c("black", "black"), legend=c("t=0", "t=1"))