今回も対象を KdV 方程式として, パラメータ $\epsilon$ をもつ次の式
$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \epsilon^2 \frac{\partial^3 u}{\partial x^3} = 0 $$を考えることにしよう．
さて、複数の保存量の保存性を同時に離散的に再現するのはなかなか難しいので、とりあえず今回は保存量のどれかに着目してその保存性を離散的に再現することを目標としよう.
前回までも示したように KdV 方程式には無限個の保存量があって簡単な方から $I_1, I_2, \cdots$ などと表されているが、今回はそのうち、
$$ I_3 := \int_0^L \left\{ \frac{u^3}{3} - \epsilon^2 (u_x)^2 \right\}dx $$に着目して、これに対応する離散量がきちんと保存されるように考えよう．
まずこのとき、$I_3 = \int_0^L G(u,u_x) dx$ として $G(u,u_x)$ を定義すると、今想定している周期的境界条件下ではその変分導関数は
$$ \frac{\delta G}{\delta u} = u^2 + 2\epsilon^2 u_{xx} $$となることが簡単な計算で確認できるので、もともとの KdV 方程式は
$$ \frac{\partial u}{\partial t} = - \frac{1}{2} \frac{\partial}{\partial x} \left( \frac{\delta G}{\delta u} \right) $$と書けることがわかる． すると、この形からほぼ自動的に $I_3$ が保存量であることが以下のようにして理解できる(境界条件は周期的境界条件としている)．
\begin{eqnarray} \frac{d}{dt} I_3[u] & = & \int_0^L \frac{\delta G}{\delta u} \frac{\partial u}{\partial t} dx \\ & = & \int_0^L \left( \frac{\delta G}{\delta u} \right) \left( -\frac{1}{2} \frac{\partial }{\partial x} \right) \left( \frac{\delta G}{\delta u} \right) dx \\ & = & - \frac{1}{2} \int_0^L \left( \frac{\delta G}{\delta u} \right) \frac{\partial }{\partial x} \left( \frac{\delta G}{\delta u} \right) dx = 0 \end{eqnarray}要は、方程式が
$$ \frac{\partial u}{\partial t} = A \left( \frac{\delta G}{\delta u} \right) $$という格好をしていて、作用素 $A$ が歪対称ならば $\int G dx$ は基本的に保存量だし、作用素 $A$ が定値性をもつならば $\int G dx$ は基本的に散逸量だ、という簡単な原理が有って、これはその一例だと言っているだけだ．
まず、このとき、保存したい保存量に対応する離散量を定義しよう．これは離散近似になっていればなんでも良いので、例えば
$$ I_{d3}[u] := {\sum_{k = 0}^N}'' G_d(u)_k \Delta x, $$ただし、 $$ G_d(u)_k := \frac{u_k^3}{3} - \epsilon^2 \left( \frac{ (\delta^{+} u_k)^2 + (\delta^{-} u_k)^2 }{2} \right) $$
としてみよう． ただし、$\sum''$ は端だけ重みを半分にする和、つまり、いわゆる台形則で、 $ \sum_{k = 0}^N '' f_k := \frac{f_0}{2} + f_1 + f_2 + \cdots + f_{N-1} + \frac{f_N}{2} $ という定義である． ただし境界条件は離散的な周期的境界条件
$$ \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. $$を想定、設定している．
次に行うのは、「離散」変分導関数の計算だ．詳細は省くが(時間があれば白板で説明する)、以下のようになるだろう．
$$ \frac{\delta G_d}{\delta (u,v)} = \frac{ u^2 + u v + v^2 }{3} + 2 \epsilon^2 \delta_k^{(2)} \left( \frac{u + v}{2} \right) $$そして、最終的に次のように差分スキームを「設計」することになる．
$$ \frac{ u^{+} - u }{ \Delta t} = - \frac{1}{2} \delta_k^{(1)} \frac{\delta G_d}{\delta (u^{+},u)} $$ただし、$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} + \frac{1}{2} \delta_k^{(1)} \left( \frac{ (u^{+})^2 + u^{+} u + u^2 }{3} \right) + \epsilon^2 \delta_k^{(3)} \left( \frac{u^{+} + u}{2} \right) = 0 $$ただし、$n = 0,1,2, \cdots$, $k = 0, 1, \cdots, N-1$ である．
では早速計算へ向けて準備を始めよう．
まずは 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
avt = ( up.^2 + up .* u + u.^2 )/6
return up - u + Δt * ( diff1(avt) + ϵ^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 というパッケージを使おう．
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 ]
保存量 $I_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 法による結果や前回の単純な差分スキームの結果とよく似ているようにも見えるが、保存量 $I_3$ が大変良く保存されている点が大きく異なる ことに注目しよう．
さらに、(たぶん「ついで」的な副作用だと思われるが) 保存量 $I_2$ の「ぶれ」もこれまでのスキームの半分程度で済んでいるし、保存量 $I_1$ はもともといずれのスキームでも保存生が良いし、ということで、いまのところ、保存量 $I_1, I_2, I_3$ に着目する限り、今回のスキームが最も「良い」ように思われる．
KdV 方程式に対して、初期値を適当に 3つのソリトンで近似した状態から計算をはじめた場合、保存量 $I_1, I_2, I_3$ の保存性が実際にどうなるか、各スキームで確認してみよ． ただしその際、ソリトンの「追い抜き」が二回以上起こるような長時間計算を行うこととせよ．