ここでもやはり、Koreweg de Vries 方程式として, パラメータ $\epsilon$ をもつ次の式
$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \epsilon^2 \frac{\partial^3 u}{\partial x^3} = 0 $$を考えることにしよう．
詳細は別資料に記載してあるが、ここでは解を区分的一次多項式(つまり折れ線グラフだ)で近似することにして、それぞれのポイントの高さを集めたベクトルを $u$ として、次のような半離散スキームを得る．
$$ \Phi \frac{d u}{dt} = - \Psi(u) u - \epsilon^2 D_1 \Phi^{-1} D_2 u $$ただし、$\Phi, \Psi(u), D_1, D_2$ はそれぞれ行列である．
この数値スキームによってベクトル $u$ に対する時間微分 $du/dt$ が得られるので、あとは Runge-Kutta 法などを使えばよいだけだ．
では早速計算してみよう．常微分方程式の解法には 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.0001 # あとでやってみるとわかるが、今回のスキームはちょっと繊細で、Δt を細かくしないとうまく動かない．
Zabusky and Kruskal https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.15.240 で使われている初期値を設定しよう．
u0 = [ cos(π*k*Δx) for k in 0:N-1 ]
念のためにプロットして確かめよう
using Plots
gr()
plot(u0)
スキームに出現する行列の類を作っていく
v = ones(N-1)
Φ = (2/3)*eye(N) + (1/6)*( diagm(v, 1) + diagm(v, -1))
Φ[1,N] = Φ[N,1] = 1/6
Φ = Δx * Φ
function cyclic_diff_center(u)
v = similar(u)
n = length(u)
v[1] = u[2] - u[n]
v[n] = u[1] - u[n-1]
for k in 2:n-1
v[k] = u[k+1] - u[k-1]
end
return v
end
function Ψ(u)
vc = (1/3) * cyclic_diff_center(u)
vf = (1/6) * diff(u)
corner = (1/6) * ( u[1] - u[end] )
M = diagm(vc, 0) + diagm(vf, 1) + diagm(vf, -1)
M[1,end] = M[end,1] = corner
return M
end
v = (1/2) * ones(N-1)
D1 = diagm( v, -1 ) - diagm( v, 1 )
D1[1,end] = 1/2
D1[end,1] = -1/2
v = ones(N-1)
D2 = 2 * eye(N) - diagm( v, 1 ) - diagm(v,-1)
D2[1,end] = D2[end,1] = -1
D2 = (1/Δx) * D2
これで、常微分方程式としての右辺を定義することが出来るようになった．
function f(u)
b1 = Φ \ (D2 * u)
b2 = - Ψ(u)*u - ϵ^2 * (D1 * b1)
v = Φ \ b2
return v
end
そして、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
試しに、初期値をプロットしてみる．
さて、これで準備は整ったので、早速、計算してみよう
using ProgressMeter
u = copy(u0)
sq_u = copy(u0)
@showprogress for i in 1:10000
u = RK(u)
if i % 1000 == 0 # 1000ステップおきにデータを格納するとして
sq_u = hcat(sq_u, u)
end
end
試しに結果の数字をちょろっと見てみる．
sq_u
plot( sq_u[:,11])
あまり変なことにはなっていなさそうなので、全体も同時にプロットしてみよう．
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) * Δx
I2(u) = sum( u.^2 ) * Δx
function I3(u)
r1 = sum( (u.^3)/3 )
r2 = sum( diff(u).^2 )/(Δ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: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 = [ OneSoliton(1, 0.5, k*Δx) + OneSoliton(0.5, 1.2, k*Δx) for k in 0:N-1 ]
試しにプロットして確認しておこう
plot(u0)
早速計算してみよう
u = copy(u0)
sq_u = copy(u0)
@showprogress for i in 1:50000 # この数字のままだと、計算時間が２０分以上かかるかもよ？
u = RK(u)
if i % 1000 == 0
sq_u = hcat(sq_u, u)
end
end
ではプロットしてみよう
plot(sq_u, legend = :no)
ちょっとわかりにくいので上から見よう．
X = 0:Δx:2-Δx # 真面目に x の範囲を書き…
T = 0:1000Δt:50000Δ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)
これもこれまで同様、二つのソリトンが衝突するあたりで数値誤差が発生している様子だ．
まあ、おおよそ上手くいっていると言えよう．
熱拡散方程式
$u_t = u_{xx}$
をディリクレ境界条件 ($u(0,t) = 0, u(L,t) = 100$) などの境界条件のもとで、有限要素法スキームを自分で構成して数値計算してみよう．