偏微分方程式: 差分法(2) 構造保存

差分法をベースに、構造保存数値解法を導入する

さて、前回は差分法によって無数のバリエーションを作れるスキームのうち、もっとも単純なものののひとつである、単なる平均的なスキームを設計して実際に動かしてみた．確かに無事に動いたので、それはそれでよいのだが、保存量の保存性があまりよくないのはちょっと気になるところであった．

そこで今回は、常微分方程式の話の “微分方程式の不変量: Volterra モデル” の回で扱ったような構造保存数値解法を導入して、偏微分方程式でもこの考え方がうまくいくか調べてみよう．

対象は KdV 方程式で

今回も対象を 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}, \cr \displaystyle u_{-1} & := & u_{N-1}, \cr \displaystyle u_{N} & := & u_{0}, \cr \displaystyle u_{N+1} & := & u_{1} \cr \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 で使われているパラメータ等にあわせよう．

1

const ϵ,L = 0.022, 2

N や他のパラメータはとりあえず次のような感じで．

1
2
3


N = 200
Δx = L/N
Δt = 0.001

配列の添字を自由にしたいので、OffsetArrays package を用いるよ．

1

using OffsetArrays

まず、境界条件を適用するシーンが何回かありそうなので、関数にしておこう．

1
2
3
4
5
6
7


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

pBC! (generic function with 1 method)

型を一回実際に使っておくと、あとで楽なので、初期値を作ってしまおう．まず、

1

u0 = OffsetArray(Float64, -2:N+1);

としてから、Zabusky and Kruskal https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.15.240 で使われている初期値を設定しよう．

1
2


u0[0:N-1] = [ cos(π*k*Δx) for k in 0:N-1 ]
pBC!(u0)

OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}} with indices -2:201: 0.998027 0.999507 1.0
0.999507 0.998027 0.995562 0.992115 0.987688 0.982287 0.975917 0.968583 0.960294 0.951057 ⋮
0.951057 0.960294 0.968583 0.975917 0.982287 0.987688 0.992115 0.995562 0.998027 0.999507 1.0
0.999507

一階から三階まで、差分作用素を関数として定義しておこう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


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$ に合わせたいので．

1
2
3
4
5


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

f (generic function with 1 method)

この式がゼロになるような解 up を「探す」わけだが、それには NLsolve package が使える．しかし、NLsolve の求解は通常の配列しか受け付けないので、通常の、外へはみ出てない配列と offsetarray を変換する関数を NLsolve 用に作っておく．

1
2
3
4
5
6
7
8


toArray(u) = u[0:N-1] # こっちは簡単．というか、定義しなくてもいいぐらい．

function toOffset(u)
    ou = similar(u0) # u0 と同じ型だしね．
    ou[0:N-1] = u #  はみ出ていない中身は本質的に同じ．
    pBC!(ou)
    return ou
end 

toOffset (generic function with 1 method)

これを使って、数値スキームを NLsolve で使える形に直しておこう．

1
2
3
4


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 (generic function with 1 method)

これで良い．この f_solve 関数がゼロになるように数値解を探せば良いことになる．さて、ではいつものように NLsolve を使う準備をしよう．

1
2
3
4
5
6
7
8
9


using NLsolve

function nls(func, params...; ini = [0.0])
    if typeof(ini) <: Number
        return nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini]).zero[1]
    else
        return nlsolve((vout,vin)->vout .= func(vin,params...), ini).zero
    end
end

nls (generic function with 1 method)

試しに 1ステップだけ解いてみよう．

1

u1 = nls(f_solve, toArray(u0), ini= toArray(u0) )

200-element Array{Float64,1}: 0.999995 0.9996
0.998218 0.99585 0.992499 0.988167 0.982858 0.976578 0.969332 0.961127 0.951971 0.941873 0.930843 ⋮
0.928705 0.939882 0.950135 0.959453 0.967827 0.975248 0.981708 0.987201 0.991722 0.995265 0.997826 0.999404

図にしてみて、変でないことをたしかめておこう．

1
2
3
4


using Plots
gr()

plot( u1 )

svg

ふむ、こんな感じか．

さて、これで準備は整ったので、早速、計算してみよう以前に説明したように、ProgressMeter というパッケージを使おう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14


using ProgressMeter

u = toArray(u0)
sq_u = copy( u )

# 1step だけ、推測値を現在の値にして, 計算 & 更新
u, u_old  = nls( f_solve, u, ini = u ), u

@showprogress for i in 2:1000
   u, u_old = nls( f_solve, u, ini = 2*u - u_old ), u  
   if i % 100 == 0  # 100ステップおきにデータを格納するとして
        sq_u = hcat(sq_u, u)
    end
end

Progress: 100%|█████████████████████████████████████████| Time: 0:02:53

試しに結果の数字をちょろっと見てみる．

1

sq_u

200×11 Array{Float64,2}: 1.0 0.956188 0.859773 … -0.523497 -0.65372 -0.66043
0.999507 0.964314 0.871821 -0.56331 -0.622243 -0.596428 0.998027 0.971693 0.883488 -0.545669 -0.543084 -0.489377 0.995562 0.978299 0.894758 -0.466434 -0.407213 -0.32827
0.992115 0.984108 0.905613 -0.316858 -0.202574 -0.100789 0.987688 0.989094 0.916032 … -0.0859584 0.0816345 0.201586 0.982287 0.993231 0.925998 0.233197 0.446155 0.573989 0.975917 0.996491 0.935487 0.631506 0.867639 0.984824 0.968583 0.998847 0.944478 1.06779 1.28511 1.3664
0.960294 1.00027 0.952947 1.45909 1.60219 1.62503
0.951057 1.00073 0.960869 … 1.69739 1.72197 1.6806
0.940881 1.00021 0.968216 1.70126 1.60281 1.5148
0.929776 0.99866 0.974959 1.46902 1.28717 1.18363
⋮ ⋱ ⋮
0.929776 0.808511 0.690693 0.414712 1.11107 1.28326
0.940881 0.823902 0.706243 0.72858 1.28025 1.11804
0.951057 0.838795 0.721563 … 1.01544 1.28168 0.83534
0.960294 0.853174 0.736643 1.2081 1.11449 0.503781 0.968583 0.867019 0.751475 1.25092 0.827434 0.182433 0.975917 0.880312 0.766046 1.12904 0.490884 -0.0935237 0.982287 0.893033 0.780347 0.877129 0.165582 -0.31119
0.987688 0.905162 0.794365 … 0.55918 -0.112196 -0.471763 0.992115 0.916679 0.808088 0.237942 -0.3289 -0.582562 0.995562 0.927561 0.821504 -0.0440674 -0.485122 -0.651883 0.998027 0.937789 0.8346 -0.266818 -0.58726 -0.686221 0.999507 0.947339 0.847361 -0.425784 -0.642087 -0.689086

あまり変なことにはなっていなさそうなので、プロットしてみよう．

1
2


# まず念のために初期値を描いてみる
plot(sq_u[:,1]) 

svg

1

plot(sq_u[:,11]) # 最後の結果を描いてみる

svg

1

plot(sq_u)  # 一応、全部重ねて描いてみる

svg

以前同様、ソリトンへの分解が見られる．

ソリトンはその「高さ」によって移動速度が異なるため(背の高いほうが速い)、時間発展に従って速いものが前に素早く移動し、遅いものはだんだんとりのこされる．いま計算した結果は、だんだん分離するソリトンを見ていることになる．

wikipedia https://ja.wikipedia.org/wiki/KdV%E6%96%B9%E7%A8%8B%E5%BC%8F の右下に、同じ計算を motion gif にしたものが掲載されているので見てみよう．

今回も KdV 方程式の保存量をチェックしよう

KdV 方程式の無限個ある保存量のうち、簡単な方からいくつか

$$ \left\{ \begin{array}{rcl} \displaystyle I_1 & := & \int_0^L u dx, \cr \displaystyle I_2 & := & \int_0^L u^2 dx, \cr \displaystyle I_3 & := & \int_0^L \left\{ u^3 / 3 - \epsilon^2 (u_x)^2 \right\}dx, \cr \displaystyle & \vdots & \end{array} \right. $$

をチェックしよう．

まずそれぞれをプログラムしておく．

1
2
3
4
5
6
7
8
9


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

早速、それぞれの量の時間発展を計算してみよう．

1

sq_I1 = [ I1(sq_u[:,k]) for k in 1:11 ]

11-element Array{Float64,1}: -2.44249e-17 -9.54792e-17 -5.39568e-16 -1.34226e-15 -6.9833e-16 -6.36158e-16 -5.62883e-16 8.31557e-16 2.59959e-15 5.2558e-15 8.44658e-15

1

sq_I2 = [ I2(sq_u[:,k]) for k in 1:11 ]

11-element Array{Float64,1}: 1.0
0.999987 0.999927 0.999596 0.997268 0.992593 0.988988 0.987213 0.986725 0.986971 0.987608

1

sq_I3 = [ I3(sq_u[:,k]) for k in 1:11 ]

11-element Array{Float64,1}: -0.0047765 -0.0047765 -0.0047765 -0.0047765 -0.0047765 -0.0047765 -0.0047765 -0.0047765 -0.0047765 -0.0047765 -0.0047765

保存量 $I_3$ が大変良く保存されている ことに着目しよう. これは狙い通りと言って良さそうだ．

複数のソリトンの挙動を追いかけてみよう

以前同様、二つのソリトンを並べた初期値を用意して、速いほうが遅い方を追い越すシチュエーションを計算してみよう．このとき、phase shift とよばれる現象(追い越すほうが前へ、追い越される方が後ろへジャンプする)が発生するとされているので、そうしたものがきちんと見えるか、注意しよう．

疑似 2-ソリトン解

高さ $H$, 中心 $x_0$ の 1-ソリトン解は次の式で書ける．

$$ u(x) = H sech^2( \sqrt{\frac{H}{12\epsilon^2}} (x-x_0) ) $$

そこで、少し離してこれを二つ配置することで 2-ソリトンを擬似的に初期配置してみよう．

1
2
3
4
5
6
7


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)

OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}} with indices -2:201: 1.035e-6
8.59703e-7 8.00653e-6 1.04091e-5 1.35327e-5 1.75936e-5 2.28731e-5 2.97369e-5 3.86603e-5 5.02614e-5 6.53436e-5 8.49514e-5 0.000110443 ⋮
4.56736e-6 3.7938e-6
3.15126e-6 2.61754e-6 2.17421e-6 1.80597e-6 1.5001e-6
1.24603e-6 1.035e-6
8.59703e-7 8.00653e-6 1.04091e-5

試しにプロットして確認しておこう

1

plot( toArray(u0) )

svg

早速計算してみよう

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


Δ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

Progress: 100%|█████████████████████████████████████████| Time: 0:07:16

ではプロットしてみよう

1

plot(sq_u, legend = :no)

svg

ちょっとわかりにくいので上から見よう．

1
2
3


X = Δx:Δx:2        # 真面目に x の範囲を書き…
T = 0:100Δt:2500Δt    # 真面目に t の範囲を書いて、
plot(X, T, sq_u')

svg

下側の辺が初期値で、全体として上へ時間発展している図になる．

やはり、二つのソリトンが基本的に一定速度で動いていること、背の高い、速いソリトンが後ろから追い越すときに、phase shift が起きていること、などがきちんと見て取れる．

この計算では保存量は？

1

sq_I1 = [ I1(sq_u[:,k]) for k in 1:26 ]

26-element Array{Float64,1}: 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198

1

sq_I2 = [ I2(sq_u[:,k]) for k in 1:26 ]

26-element Array{Float64,1}: 0.137543 0.137543 0.137544 0.137544 0.137545 0.137545 0.137546 0.137547 0.13755 0.137554 0.137562 0.137575 0.137595 0.137625 0.137662 0.137698 0.137717 0.137708 0.137677 0.137639 0.137605 0.137581 0.137565 0.137556 0.13755 0.137548

1

sq_I3 = [ I3(sq_u[:,k]) for k in 1:26 ]

26-element Array{Float64,1}: 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467 0.0239467

1

plot(sq_I1 .- 0.260198 , marker = :circle) # おおよその値を引いてズレをみる

svg

1

plot( sq_I2 .- 0.137543, marker = :circle)

svg

1

plot( sq_I3 .- 0.0239467, marker = :circle)

svg

二つのソリトンが衝突するあたりで数値誤差が発生している様子など、おおざっぱな意味では前々回の method of line + Runge-Kutta 法による結果や前回の単純な差分スキームの結果とよく似ているようにも見えるが、保存量 $I_3$ が大変良く保存されている点が大きく異なる ことに注目しよう．

さらに、(たぶん「ついで」的な副作用だと思われるが) 保存量 $I_2$ の「ぶれ」もこれまでのスキームの半分程度で済んでいるし、保存量 $I_1$ はもともといずれのスキームでも保存生が良いし、ということで、いまのところ、保存量 $I_1, I_2, I_3$ に着目する限り、今回のスキームが最も「良い」ように思われる．

Report No.17 疑似 3-ソリトン

KdV 方程式に対して、初期値を適当に 3つのソリトンで近似した状態から計算をはじめた場合、保存量 $I_1, I_2, I_3$ の保存性が実際にどうなるか、各スキームで確認してみよ．ただしその際、ソリトンの「追い抜き」が二回以上起こるような長時間計算を行うこととせよ．

Report No.18 なぜうまくいくのか

この数値スキームだとなぜ保存量 $I_3$ がこんなにも良く保存されるのだろうか．考えてみよう．