前回紹介した KdV 方程式として, パラメータ $\epsilon$ をもつ次の式
$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \epsilon^2 \frac{\partial^3 u}{\partial x^3} = 0 $$を考えることにしよう．
差分法は微分を差分、つまり、Taylor 展開で近似するものである． 近似できていれば良いのだ！ と考えれば、こんなに自由な方法もないわけで、まあやりたい放題といって良い．
そこで、たとえば今回は、あちこち平均したり差分したりして作った以下の離散式(数値スキーム)、をもとの KdV 方程式の近似式として考えよう． ただし、$u_k^{(n)} \cong u(k\Delta x, n\Delta t)$ として、下の式は簡潔に書くために $u = u_k^{(n)}$, $u^{+} = u_k^{(n+1)}$ と記している．
$$ \frac{ u^{+} - u }{\Delta t} + \left( \frac{u^{+} + u}{2} \right) \delta_k^{(1)} \left( \frac{u^{+} + u}{2} \right) + \epsilon^2 \delta_k^{(3)} \left( \frac{u^{+} + u}{2} \right) = 0 $$ただし、$n = 0,1,2, \cdots$, $k = 0, 1, \cdots, N-1$ で、離散境界条件は以前と同じように周期的境界条件を離散化した
$$ \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. $$としよう．
では早速計算へ向けて準備を始めよう．
まずは 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
まず、境界条件を適用するシーンが何回かありそうなので、関数にしておこう．
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
型を一回実際に使っておくと、あとで楽なので、初期値を作ってしまおう．まず、
u0 = OffsetArray(Float64, -2:N+1);
としてから、Zabusky and Kruskal https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.15.240 で使われている初期値を設定しよう．
u0[0:N-1] = [ cos(π*k*Δx) for k in 0:N-1 ] # 中身は普通に作って、
pBC!(u0) # 境界条件を満たすようにする
一階から三階まで、差分作用素を関数として定義しておこう．
function diff1(u)
v = similar(u)
v[0:N-1] = [( u[k+1] - u[k-1] )/(2Δx) for k in 0:N-1 ]
return pBC!(v)
end
function diff2(u)
v = similar(u)
v[0:N-1] = [ ( u[k+1] -2u[k] + u[k-1] ) /(Δx^2) for k in 0:N-1 ]
return pBC!(v)
end
diff3(u) = diff1(diff2(u))
ただし、全体に $\Delta t$ をかけておく．式全体のオーダーを $u$ に合わせたいので．
function f(up, u)
av = (up+u)/2
return up - u + Δt * ( av .* diff1(av) + ϵ^2 * diff3(av) )
end
この式がゼロになるような解 up を「探す」わけだが、それには NLsolve package が使える． しかし、NLsolve の求解は通常の配列しか受け付けないので、 通常の、外へはみ出てない配列と offsetarray を変換する関数を NLsolve 用に作っておく．
toArray(u) = u[0:N-1] # こっちは簡単．というか、定義しなくてもいいぐらい．
function toOffset(u)
ou = similar(u0) # u0 と同じ型だしね．
ou[0:N-1] = u # はみ出ていない中身は本質的に同じ．
pBC!(ou)
return ou
end
これを使って、数値スキームを NLsolve で使える形に直しておこう．
function f_solve(up_in, u_in) # up_in も u_in も添字が 1:N-1 の普通の配列だとしよう．
up, u = toOffset(up_in), toOffset(u_in)
return toArray( f(up,u) ) # 出力も普通の配列にして…
end
これで良い．この f_solve 関数がゼロになるように数値解を探せば良いことになる． さて、では、NLsolve を使う準備をしよう．
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
return nlsolve((vin,vout)->vout[1]=func(vin[1],params...), [ini]).zero[1]
else
return nlsolve((vin,vout)->vout .= func(vin,params...), ini).zero
end
end
試しに 1ステップだけ解いてみよう．
u1 = nls(f_solve, toArray(u0), ini= toArray(u0) )
図にしてみて、変でないことをたしかめておこう．
using Plots
gr()
plot( u1 )
ふむ、こんな感じか．
さて、これで準備は整ったので、早速、計算してみよう
このスキームは非線形方程式を解くので、時間がかかることが想定される． そこで、計算がどれくらい進んでいるのか、リアルタイムで見えるようにしよう．それには、ProgressMeter というパッケージを使えば良い(インストールしていない人は、Pkg.add コマンドでインストールしておこう)．
これを使うと、@showpprogress 時間間隔(秒) メッセージ などを for 文の頭につけるだけで、その for 文の実行状況と、残り推定時間を指定した時間間隔毎に更新して教えてくれる．
using ProgressMeter
u = toArray(u0) # メタ的視点から見ると、offsetarray は見えない格好になった．よって、ここは通常配列で．
sq_u = copy( u )
# 1step だけ、推測値を現在の値にして, 計算 & 更新
u, u_old = nls( f_solve, u, ini = u ), u
@showprogress for i in 2:1000 # 環境にもよるが、10分ぐらいかかるか？
# 授業時は 2:200 ぐらいにしておいた方が良さ気．
u, u_old = nls( f_solve, u, ini = 2*u - u_old ), u
# 2step 目から、推測値を「過去と今から計算」して入力．こうすると何割か早い．
if i % 100 == 0 # 100ステップおきにデータを格納するとして
sq_u = hcat(sq_u, u)
end
end
試しに結果の数字をちょろっと見てみる．
sq_u
あまり変なことにはなっていなさそうなので、プロットしてみよう．
# まず念のために初期値を描いてみる
plot(sq_u[:,1])
plot(sq_u[:,11]) # 最後の結果を描いてみる
plot(sq_u) # 一応、全部重ねて描いてみる
ソリトンはその「高さ」によって移動速度が異なるため(背の高いほうが速い)、時間発展に従って速いものが前に素早く移動し、遅いものはだんだんとりのこされる．いま計算した結果は、だんだん分離するソリトンを見ていることになる．
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) * Δx
I2(u) = sum( u.^2 ) * Δx
function I3(u)
r1 = sum( (u.^3)/3 )
r2 = sum( (vcat(diff(u), u[1] - u[end] )).^2 )/(Δx^2)
return (r1 - ϵ^2 * r2)*Δx
end
早速、それぞれの量の時間発展を計算してみよう．
sq_I1 = [ I1(sq_u[:,k]) for k in 1:11 ]
sq_I2 = [ I2(sq_u[:,k]) for k in 1:11 ]
sq_I3 = [ I3(sq_u[:,k]) for k in 1:11 ]
とりあえず、この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
u0[0:N-1] = [ OneSoliton(1, 0.5, k*Δx) + OneSoliton(0.5, 1.2, k*Δx) for k in 0:N-1]
pBC!(u0)
試しにプロットして確認しておこう
plot( toArray(u0) )
早速計算してみよう
Δt = 0.002 # 計算時間がかかるので、時間ステップ間隔を大きくしてしまえ．
u = toArray(u0) # メタ的視点から見ると、offsetarray は見えない格好になった．よって、ここは通常配列で．
sq_u = copy( u )
# 1step だけ、推測値を現在の値にして, 計算 & 更新
u, u_old = nls( f_solve, u, ini = u ), u
@showprogress for i in 2:2500 # 環境にもよるが、15分ぐらいかかるか？
# 授業時は 2:200 ぐらいにしておいた方が良さ気．
u, u_old = nls( f_solve, u, ini = 2*u - u_old ), u # 2step 目から、推測値を「過去と今から計算」して入力．
if i % 100 == 0 # 100ステップおきにデータを格納するとして
sq_u = hcat(sq_u, u)
end
end
ではプロットしてみよう
plot(sq_u, legend = :no)
ちょっとわかりにくいので上から見よう．
X = Δx:Δx:2 # 真面目に x の範囲を書き…
T = 0:100Δt:2500Δt # 真面目に t の範囲を書いて、
plot(X, T, sq_u')
下側の辺が初期値で、全体として上へ時間発展している図になる．
やはり、二つのソリトンが基本的に一定速度で動いていること、 背の高い、速いソリトンが後ろから追い越すときに、phase shift が起きていること、などがきちんと見て取れる．
sq_I1 = [ I1(sq_u[:,k]) for k in 1:26 ]
sq_I2 = [ I2(sq_u[:,k]) for k in 1:26 ]
sq_I3 = [ I3(sq_u[:,k]) for k in 1:26 ]
plot(sq_I1 .- 0.260198 , marker = :circle) # 変化が小さいので、おおよその値を引いて、そこからのズレをみる
plot( sq_I2 .- 0.137543, marker = :circle)
plot( sq_I3 .- 0.0239467, marker = :circle)
二つのソリトンが衝突するあたりで数値誤差が発生している様子など、先週の method of line + Runge-Kutta 法による結果とよく似ている．
KdV 方程式に対して、差分法によって他の数値スキームを作って、実際に計算して保存量などをチェックしてみよう． 特に、「保存量がなるべく保たれるような数値スキーム」を作れないか、工夫してみよう．