例えば、独立変数が空間 $x$ と時間 $t$ の時間発展問題を記述する偏微分方程式は一般に
$$ F(u, \partial u/\partial t, \partial u/\partial x, \cdots) = 0 $$と書けるが、これを大変多くの式が連立する常微分方程式
$$ \frac{\partial u}{\partial t} = f(u, V, W, \cdots) $$に帰着ないしは近似させて、この連立常微分方程式を数値計算する方法がいわゆる "method of line" である． (線の方法などと直訳されるが、口語っぽすぎて訳としてわかりにくいかも)
これだけだとよくわからないだろうから、もう少し例を交えて解説すると、たとえば、 $ u(x,t) $ という関数を, $x$ 方向だけ $ u_k(t) \cong u(k\Delta x, t) $ と離散近似化して、 $ u(x,t) $ 全体の代わりに(近似で) $ \{ u_0(t), u_1(t), u_2(t), \cdots \}$ を扱うことにするのだ．
そして、例えば $\partial u/\partial x$ を差分近似 $ ( u_{k+1} - u_{k-1}) / 2\Delta x$ などで置き換えるなどの処理を行えば、上の $f$ を作ることが出来る．
そして、できあがった常微分方程式の数値計算には、これまで学んだ Runge-Kutta 法や線形多段解法、その他もろもろ、どれでも好きな方法を使えば良い．
Koreweg de Vries 方程式、略して KdV 方程式は、非線形な偏微分方程式であるにも関わらず、進行する孤立波を解として持つという、大変興味深い特徴を持つ方程式だ． ここでは、KdV 方程式として, パラメータ $\epsilon$ をもつ次の式
$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \epsilon^2 \frac{\partial^3 u}{\partial x^3} = 0 $$を考えることにしよう．
ここでは、空間領域 $\Omega = [0, L]$ に対し、その分割数を $N$ とし、境界条件として周期的境界条件を設定しよう．
よって、$\Delta x = L/N$ として従属変数を $ U(t) = \{ u_0(t), u_1(t), \cdots , u_{N-1}(t)\}$ とし、境界条件は たとえば仮に外側にはみ出す点を仮想的に考えたとして $$ \left\{ \begin{array}{rcl} \displaystyle u_{-2} & := & u_{N-2}, \\ \displaystyle u_{-1} & := & u_{N-1}, \\ \displaystyle u_{N} & := & u_{0}, \\ \displaystyle u_{N+1} & := & u_{1} \\ \end{array} \right. $$ と解釈することになる．
こうしておいて、空間微分をたとえば差分近似で置き換える． 一階微分の近似である一階中心差分ならば
$$ \delta_k^{(1)} u_k := \frac{ u_{k+1} - u_{k-1} }{2\Delta x} $$だし、二階中心差分ならば
$$ \delta_k^{(2)} u_k := \frac{ u_{k+1} - 2u_k + u_{k-1} }{\Delta x^2} $$だし、さらに三階などは
$$ \delta_k^{(3)} := \delta_k^{(1)}\delta_k^{(2)} $$となる． これらを使った場合は、全体は
$$ \frac{d u_k(t)}{dt} = - u_k (t) \delta_k^{(1)} u_k(t) - \epsilon^2 \delta_k^{(3)} u_k(t), k = 0,1,\cdots,N-1 $$という N元連立常微分方程式で近似できることになる．
では早速計算してみよう．常微分方程式の解法には Runge-Kutta 法を使うことにしよう．
まずは wikipedia https://ja.wikipedia.org/wiki/KdV%E6%96%B9%E7%A8%8B%E5%BC%8F でも引用されている Zabusky and Kruskal の1965年の古い論文 https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.15.240 で使われているパラメータ等にあわせよう．
const ϵ,L = 0.022, 2
N や他のパラメータはとりあえず次のような感じで．
N = 200
Δx = L/N
Δt = 0.001
(インストールしてない人はインストールしておこう)
using OffsetArrays
すると次のような感じで、自由に決めた添字で配列を使える．
u0 = OffsetArray(Float64, -2:N+1);
境界条件を適用するシーンが何回かありそうなので、関数にしておこう．
function pBC!(u) # 境界条件に従わせる.
u[-1] = u[N-1]
u[-2] = u[N-2]
u[N] = u[0]
u[N+1] = u[1]
return u
end
一階から三階まで、差分作用素を関数として定義しておこう．
function diff1(u)
v = copy(u)
for k in 0:N-1
v[k] = ( u[k+1] - u[k-1] )/(2Δx)
end
return pBC!(v)
end
function diff2(u)
v = copy(u)
for k in 0:N-1
v[k] = ( u[k+1] -2u[k] + u[k-1] ) /(Δx^2)
end
return pBC!(v)
end
diff3(u) = diff1(diff2(u))
これで、常微分方程式としての右辺を定義することが出来るようになった．
f(u) = - u .* diff1(u) - ϵ^2 * diff3(u) # .* で、配列の要素ごとの掛け算が出来る．
そして、Runge-Kutta 法を定義する．
function RK(u)
f1 = f(u)
f2 = f(u + Δt/2 * f1)
f3 = f(u + Δt/2 * f2)
f4 = f(u + Δt * f3)
return u + Δt * (f1 + 2*f2 + 2*f3 + f4)/6
end
Zabusky and Kruskal https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.15.240 で使われている初期値を設定しよう．
for k in 0:N-1
u0[k] = cos(π*k*Δx)
end
pBC!(u0)
試しに、初期値をプロットしてみる．
using Plots
gr()
plot(u0[1:N]) # plot は普通の配列しか理解できないので、そこだけ切り取って渡す
ふむ、こんな感じか．
さて、これで準備は整ったので、早速、計算してみよう
u = copy(u0)
sq_u = copy(u0[1:N]) # 記録しておくのは、内点の情報だけでいい．通常の配列扱いできるし．
for i in 1:1000
u = RK(u)
if i % 200 == 0 # 200ステップおきにデータを格納するとして
sq_u = hcat(sq_u, u[1:N])
end
end
試しに結果の数字をちょろっと見てみる．
sq_u
あまり変なことにはなっていなさそうなので、プロットしてみよう．
# まず念のために初期値を描いてみる
plot(sq_u[:,1])
plot(sq_u[:,6]) # 最後の結果を描いてみる
plot(sq_u) # 一応、全部重ねて描いてみる
実は、KdV 方程式の解は複数のソリトンの重ね合わせ的なもので書けることがわかっている．
そして、ソリトンはその「高さ」によって移動速度が異なるため(背の高いほうが速い)、時間発展に従って速いものが前に素早く移動し、遅いものはだんだんとりのこされる．いま計算した結果は、だんだん分離するソリトンを見ていることになる．
wikipedia https://ja.wikipedia.org/wiki/KdV%E6%96%B9%E7%A8%8B%E5%BC%8F の右下に、同じ計算を motion gif にしたものが掲載されているので見てみよう．
上の計算結果を、どうにかして自分で動画にしてみよう．
KdV 方程式には、無限個の保存量があることがわかっている．単純な方から挙げてみると、
$$ \left\{ \begin{array}{rcl} \displaystyle I_1 & := & \int_0^L u dx, \\ \displaystyle I_2 & := & \int_0^L u^2 dx, \\ \displaystyle I_3 & := & \int_0^L \left\{ u^3/3 - \epsilon^2 (u_x)^2 \right\}dx, \\ \displaystyle & \vdots & \end{array} \right. $$という感じだ．
それぞれ、プログラムしてみよう．
I1(u) = sum(u[1:N]) * Δx
I2(u) = sum( u[1:N].^2 ) * Δx
function I3(u)
r1 = sum( (u.^3)/3 )
r2 = sum( (u[k+1] - u[k])^2 for k in 1:N-1 )/(Δx^2)
r2 += (u[1]-u[N])^2/(Δx^2) # 面倒だけど、こいつもいれておかないとね…
return (r1 - ϵ^2 * r2)*Δx
end
早速、それぞれの量の時間発展を計算してみよう．
sq_I1 = [ I1(sq_u[:,k]) for k in 1:6 ]
sq_I2 = [ I2(sq_u[:,k]) for k in 1:6 ]
sq_I3 = [ I3(sq_u[:,k]) for k in 1:6 ]
とりあえず、この3つの保存量はだいたい保存されていることがみてとれる．
先にも書いたように、ソリトンは高いほうが速い．
そこで、二つのソリトンを並べた初期値を用意して、速いほうが遅い方を追い越すシチュエーションを計算してみよう． このとき、phase shift とよばれる現象(追い越すほうが前へ、追い越される方が後ろへジャンプする)が発生するとされているので、それを観測してみよう．
高さ $H$, 中心 $x_0$ の 1-ソリトン解は次の式で書ける．
$$ u(x) = H sech^2( \sqrt{\frac{H}{12\epsilon^2}} (x-x_0) ) $$そこで、少し離してこれを二つ配置することで 2-ソリトンを擬似的に初期配置してみよう．
function OneSoliton(height, center, x) # 1-ソリトン解
d = sqrt(height/12) / ϵ
return height * sech(d*(x-center))^2
end
for k in 0:N-1
u0[k] = OneSoliton(1, 0.5, k*Δx) + OneSoliton(0.5, 1.2, k*Δx)
end
pBC!(u0)
試しにプロットして確認しておこう
plot(u0[1:N])
早速計算してみよう
u = copy(u0)
sq_u = copy(u0[1:N])
for i in 1:5000
u = RK(u)
if i % 100 == 0
sq_u = hcat(sq_u, u[1:N])
end
end
ではプロットしてみよう
plot(sq_u, legend = :no)
ちょっとわかりにくいので上から見よう．
X = Δx:Δx:2 # 真面目に x の範囲を書き…
T = 0:100Δt:5000Δt # 真面目に t の範囲を書いて、
plot(X, T, sq_u')
下側の辺が初期値で、全体として上へ時間発展している図になる．
二つのソリトンが基本的に一定速度で動いていること、 背の高い、速いソリトンが後ろから追い越すときに、phase shift が起きていること、などがきちんと見て取れる．
sq_I1 = [ I1(sq_u[:,k]) for k in 1:51 ]
sq_I2 = [ I2(sq_u[:,k]) for k in 1:51 ]
sq_I3 = [ I3(sq_u[:,k]) for k in 1:51 ]
plot(sq_I1 .- 0.260198 , marker = :circle) # 変化が小さいので、おおよその値を引いて、そこからのズレをみる
plot( sq_I2 .- 0.137543, marker = :circle)
plot( sq_I3 .- 0.0239467, marker = :circle)
どうやら、二つのソリトンが衝突するあたりで数値誤差が発生している様子だ．
KdV 方程式の数値計算を、method of line + 線形多段解法(+ 予測子修正式法) でやってみよう．