偏微分方程式: method of line

偏微分方程式の数値計算, 一番単純な考え方: method of line

さて、ここからはこれまで扱ってきた常微分方程式ではなく、偏微分方程式を扱おう.

例えば、独立変数が空間 $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 \in [0,L], t \in (0, T) $ という関数を, $x$ 方向だけ $x_k := k\Delta x$, $k \in 0,1, \cdots, N$ などと離散化して 解を $u_k(t) \cong u(x_k, t) $ と離散近似化し、 $u(x,t)$ の偏微分方程式の代わりに $\{ u_0(t), u_1(t), u_2(t), \cdots, u_N(t) \}$ の常微分方程式を扱うことにするのだ.

この際、例えば $x$ 方向の微分 $\partial u/\partial x$ などはを差分近似 $ ( u_{k+1} - u_{k-1}) / 2\Delta x$ などで置き換えるなどの処理を行えば良い(もちろんいろいろな方法がある).

そして、できあがった常微分方程式の数値計算には、これまで学んだ Runge-Kutta 法や線形多段解法、その他もろもろ、どれでも好きな方法を使えば良い.

KdV 方程式を試してみよう

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 \]

を考えることにしよう.

method of line を適用する

ここでは、空間領域 $\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}, \cr \displaystyle u_{-1} & := & u_{N-1}, \cr \displaystyle u_{N} & := & u_{0}, \cr \displaystyle u_{N+1} & := & u_{1} \cr \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 で使われているパラメータ等にあわせよう.

1
const ϵ,L = 0.022, 2

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

1
2
3
N = 200
Δx = L/N
Δt = 0.001

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

Julia の配列は基本的に添え字が 1から始まるが、OffsetArrays というパッケージ(JuliaBox にはデフォルトでインストールされている)を用いると添字範囲を自由に変えられる. これを使ってプログラムを容易にしよう.

まず使用宣言.

1
using OffsetArrays

すると次のような感じで、自由に決めた添字で配列を使える.

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

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

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

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
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))

これで、常微分方程式としての右辺を定義することが出来るようになった.

1
f(u) = - u .* diff1(u) - ϵ^2 * diff3(u)  # .* で、配列の要素ごとの掛け算が出来る.

そして、いつものごとく Runge-Kutta 法を定義する.

1
2
3
4
5
6
7
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 で使われている初期値を設定しよう.

1
2
3
4
5
u0 = OffsetArray(Float64, -2:N+1);
for k in 0:N-1
    u0[k] = cos(π*k*Δx) # cos の中身は、円周率パイ * k * デルタx 
end
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
using Plots
gr()
1
plot(u0[1:N]) # plot は普通の配列しか理解できないので、そこだけ切り取って渡す

svg

ふむ、こんな感じか.

さて、これで準備は整ったので、早速、計算してみよう

1
2
3
4
5
6
7
8
9
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

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

1
sq_u

200×6 Array{Float64,2}: 0.999507 0.871751 0.696739 0.403807 -0.576518 -0.602904 0.998027 0.883419 0.712713 0.216225 -0.559204 -0.49782 0.995562 0.894689 0.724192 0.0685818 -0.479925 -0.336876 0.992115 0.905544 0.730027 -0.0236337 -0.332198 -0.110168 0.987688 0.915965 0.730365 -0.0459073 -0.100761 0.196353 0.982287 0.925932 0.726934 0.00891994 0.221211 0.577461 0.975917 0.935423 0.722938 0.150279 0.63133 1.00675 0.968583 0.944416 0.722523 0.377671 1.08832 1.41052 0.960294 0.952887 0.729886 0.681382 1.50779 1.68633 0.951057 0.96081 0.748191 1.02232 1.76394 1.73734 0.940881 0.96816 0.77846 1.33304 1.76168 1.54701 0.929776 0.974907 0.818688 1.52149 1.49928 1.18419 0.917755 0.981019 0.863473 1.51922 1.07818 0.759683 ⋮ ⋮
0.940881 0.706179 0.53463 0.152727 0.740955 1.14066 0.951057 0.721498 0.54771 0.162266 1.04428 0.835579 0.960294 0.736578 0.559752 0.216979 1.25063 0.489555 0.968583 0.751409 0.570674 0.315473 1.2935 0.160084 0.975917 0.765979 0.580689 0.453949 1.15864 -0.115537 0.982287 0.780279 0.590331 0.614304 0.886094 -0.331509 0.987688 0.794296 0.600393 0.77233 0.551648 -0.487252 0.992115 0.808019 0.611766 0.891584 0.219991 -0.595461 0.995562 0.821435 0.625202 0.94187 -0.0644904 -0.660933 0.998027 0.83453 0.641055 0.902916 -0.287029 -0.694332 0.999507 0.847291 0.659049 0.782059 -0.442896 -0.695195 1.0 0.859704 0.67817 0.602996 -0.538929 -0.667305

あまり変なことにはなっていなさそうなので、プロットしてみよう. まず念のために初期値を描いてみる

1
plot(sq_u[:,1]) 

svg

最後の結果を描いてみる

1
plot(sq_u[:,6])

svg

一応、全部重ねて描いてみる

1
plot(sq_u)

svg

何が起きているか?

実は、KdV 方程式の解は複数のソリトン(おおよそ、伝搬していく孤立波のことと思えば良い)の重ね合わせ的なもので書けることがわかっている.

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

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
10
11
12
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

プログラムを修正して数値計算をやり直してもよいが、すでに計算結果はあるのでそれをもとにそれぞれの量の時間発展を計算してみよう.

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

6-element Array{Float64,1}: -3.10862e-17 -1.33227e-17 -8.54872e-17 -2.19824e-16 -2.5091e-16 -2.81997e-16

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

6-element Array{Float64,1}: 1.0
1.00015 1.00576 1.02334 1.02791 1.02588

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

6-element Array{Float64,1}: -0.0047765
-0.00477863 -0.00395926 -0.000880186 -0.00083225 -0.0013602

この数値計算だと、$I_1$ はかなりよく保存されているが、他の 2つのつの保存量は「かなり甘く」保存されていることがみてとれる.

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

先にも書いたように、ソリトンは高いほうが速い.

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

(疑似) 2-ソリトン解

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

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

1
2
3
4
function OneSoliton(height, center, x)
    d = sqrt(height/12) / ϵ
    return height * sech(d*(x-center))^2 
end

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

注: 2-ソリトン解は 1-ソリトン解を2つ並べたものとは数式も形状も異なるが、ソリトン同士がある程度離れている場合は 1-ソリトン解 2つでかなりよく近似できる. そして、2-ソリトン解の数式は面倒だ…

1
2
3
4
for k in 0:N-1
    u0[k] = OneSoliton(1, 0.5, k*Δx) + OneSoliton(0.5, 1.2, k*Δx)
end
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(u0[1:N])

svg

早速計算してみよう

1
2
3
4
5
6
7
8
9
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

ではプロットしてみよう

1
plot(sq_u, legend = :no)

svg

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

1
2
3
X = Δx:Δx:2
T = 0:100Δt:5000Δt
plot(X, T, sq_u')

svg

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

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

この計算では保存量は?

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

51-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

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

51-element Array{Float64,1}: 0.137543 0.137546 0.137547 0.137548 0.137548 0.137548 0.137548 0.137548 0.137547 0.137546 0.137546 0.137545 0.137543 ⋮
0.137389 0.13742 0.137446 0.137468 0.137486 0.137502 0.137513 0.13752 0.137527 0.137532 0.137537 0.137539

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

51-element Array{Float64,1}: 0.0239467 0.0239473 0.0239477 0.0239479 0.0239478 0.0239477 0.0239479 0.0239479 0.0239477 0.0239473 0.023947 0.0239466 0.023946 ⋮
0.0238668 0.0238831 0.0238969 0.0239082 0.0239175 0.0239251 0.0239305 0.0239345 0.0239378 0.0239403 0.0239425 0.0239437

それぞれプロットしてみよう. ただし、あまり変化はないようなので、適度に初期値を引いてみよう.

1
plot(sq_I1 - sq_I1[1], marker = :circle)

svg

WARNING: No strict ticks found WARNING: No strict ticks found WARNING: No strict ticks found

1
plot( sq_I2 - sq_I2[1], marker = :circle)

svg

1
plot( sq_I3 - sq_I3[1], marker = :circle)

svg

どうやら、二つのソリトンが衝突するあたりで数値誤差が発生しているようだ、ということがわかる.

Report No.09 KdV 方程式

  1. KdV 方程式の近似解の時間発展の様子をどうにかして自分で動画にしてみよう.
  2. KdV 方程式の数値計算を、今使った Runge-Kutta スキームの部分を 線形多段解法(+ 予測子修正式法) に替えてやってみよう.