09. Method of Line for 偏微分方程式

grid Photo by Samuel Zeller on Unsplash

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

さて、ここからはこれまで扱ってきた常微分方程式ではなく、偏微分方程式を扱おう． 偏微分方程式になると独立変数が複数になるので，問題は格段にややこしくなる．

まずは独立変数が空間 $x$ と時間 $t$ で従属変数 $u$ が時間経過で変化する時間発展問題を記述する偏微分方程式を扱うことにしよう． なにかと理解しやすいし，応用範囲も広いしね．

さて，このタイプの方程式は一般に

\begin{equation} F\left(u, \frac{\partial u}{\partial t}, \, \frac{\partial u}{\partial x}, \, \cdots \right) = 0 \end{equation}

と書けると思っていいだろう． このような問題に対して、近似操作で従属変数を見かけ上増やすなどの操作によって「時間以外」の独立変数を見かけ上消して 次式のような形，

\begin{equation} \frac{\partial u}{\partial t} = f(u, V, W, \cdots) \end{equation}

すなわち，時間に対する常微分方程式に問題を変形してから 数値計算を行う方法が, 今回学ぶ，いわゆる "method of line" である．

これだけだとよくわからないだろうから、もう少し例を交えて解説しよう．

たとえば、$x, t$ という空間変数，時間変数を独立変数にもつ関数 $ u(x,t) $ $ x \in [0,L], t \in (0, T) $ が問題の解であるとしよう．

このとき，$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) \}$ の常微分方程式を扱うことにする，という算段が "method of line" だ．

この際、例えば $x$ 方向の微分 $\partial u/\partial x$ 前進差分 $ ( u_{k+1} - u_k) / \Delta x$ で離散近似するもよし， 中心差分近似 $ ( u_{k+1} - u_{k-1}) / 2\Delta x$ で離散近似するもよし，ということで離散近似化にはいろいろな方法があるので，一口で "method of line" と言っても具体的な手法はいろいろだ．

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

KdV 方程式を試してみよう

では早速，具体的な偏微分方程式を対象にこの "method of line" を試してみよう．

ここで扱う偏微分方程式として， Koreweg de Vries 方程式をとりあげよう．

この方程式(略して KdV 方程式)は、非線形な偏微分方程式であるにも関わらず、進行する孤立波を解として持つという、大変興味深い特徴を持つ方程式だ． 具体的な表式にはすこしバリエーションがあるが， ここでは、KdV 方程式として次の式

\begin{equation} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \epsilon^2 \frac{\partial^3 u}{\partial x^3} = 0, \,\,\,\,\, t > 0, \,\,\, x \in \Omega = [0, L], \end{equation}

を考えることにしよう(ただし $x,t$ はそれぞれ 1次元値，$\epsilon$ はパラメータ)．

ただし，偏微分方程式を考察するには境界条件が必要だ(境界条件なしでは解が定まらない).

ここでは，時間方向としては，空間方向の境界条件を満たし適度に滑らかな初期値: \begin{equation} u(x, 0) = u_0(x) \end{equation} が与えられるとし，空間方向については， 境界条件として以下のような周期的境界条件 \begin{equation} u(x + mL, t) = u(x, t), \,\,\, x \in \Omega, \,\, m:\mbox{ 任意の整数} \label{eq:bc} \end{equation} を設定しよう．

method of line を適用する

さっそく method of line を適用してみよう． まずは空間変数のからむ部分の離散近似化からだね．

そうだな，空間領域 $\Omega$ を均等に分割するとしよう． 分割数を $N$ とし $\Delta x = L/N$ とすると 空間の分点を $x_k = k\Delta x$ (ただし $k$ は整数)とすることができる．

そして分点 $x_k$ の上での $u(x_k,t)$ の近似値を $u_k(t)$ と書いて， 時間 $t$ における計算対象とすべき値全体 $\boldsymbol{U}(t) = \{ u_0(t), u_1(t), \cdots , u_{N-1}(t) \}$ を考える1．

次に境界条件を離散化しよう． 境界条件の離散化には何通りかのやりかたがあるが， ここではいろいろと扱いやすい「計算対象分点の外側に仮想空間分点を用意して，境界条件によってその値を決める」 というやり方を採用してみよう．

今回は，問題の中に空間微分としては最大で 3階のものがある． これを差分法で離散近似すると分点が少なくとも 4点必要で，上手に中心差分を使った場合，左右に 2点ずつの分点が必要だ． よって，中心差分を採用することにして，仮想点として，計算対象の $\boldsymbol{U}$ が用いる分点の左側の外に 2点，右側の外に 2点必要ということになる．

よってその分点 4点における近似値 $u_{-2}$, $u_{-1}$, $u_{N}$, $u_{N+1}$ を，境界条件 $(\ref{eq:bc})$ を満たす値 \begin{equation} \left\{ \begin{array}{rcl} \displaystyle u_{-2}(t) & := & u_{N-2}(t), \cr \displaystyle u_{-1}(t) & := & u_{N-1}(t), \cr \displaystyle u_{N}(t) & := & u_{0}(t), \cr \displaystyle u_{N+1}(t) & := & u_{1}(t) \cr \end{array} \right. \end{equation} に設定して，この値を計算で参照すれば良い． これで，境界条件を離散化したことになる．

こうしておいて、空間微分をたとえば差分近似で置き換える． 一階微分の近似である一階中心差分ならば

\begin{equation} \delta_k^{(1)} u_k := \frac{ u_{k+1} - u_{k-1} }{2\Delta x} \end{equation}

という近似式になるし、二階中心差分ならば

\begin{equation} \delta_k^{(2)} u_k := \frac{ u_{k+1} - 2u_k + u_{k-1} }{\Delta x^2} \end{equation}

という形だ． そして三階微分は三階の中心差分

\begin{equation} \delta_k^{(3)} := \delta_k^{(1)}\delta_k^{(2)} \end{equation}

で近似することになるだろう．

よって，空間微分に対してこれらの差分を素直に使って作った離散近似式は，上の境界条件のもとで

\begin{equation} \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 \end{equation}

という N元連立常微分方程式になるね．

早速計算してみよう！

では早速計算してみよう．常微分方程式の解法には Runge-Kutta 法を使うことにしよう．

まずは wikipedia でも引用されている Zabusky and Kruskal の1965年の古い論文 で使われているパラメータ等にあわせよう．

1const ϵ,L = 0.022, 2

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

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

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

Julia の配列は基本的に添え字が 1から始まるが、OffsetArrays というパッケージを用いると自由に設定した添字範囲で配列にアクセスできる． これを使ってプログラムを容易にしよう．

OffsetArrays パッケージのインストール (VIORAS 環境では不要！)

自分の PC に自前で julia をインストールしているような場合はこの package を自分でインストールする必要がある. 一回だけやれば良い． 下記のようにしてインストールしよう．

  VIORAS環境ではこのパッケージのインストールは不要なので，下記のインストール作業を実行しないように！

1using Pkg
2Pkg.add("OffsetArrays")



さて，このパッケージが使えるようになっているならば，まずはいつものように使用宣言だ．

1using OffsetArrays

そしてこの OffsetArrays の使い方を説明しておこう．

基本的には，まず普通の配列を作っておいて，それに OffsetArray 命令で別の名前をつけるのだ． その別名をつけるときにアクセスする添字範囲を異なるものに設定できる，という仕組みだ．

だから，名前は複数あっても実体は 1つしかない． もとの名前の配列の中身を書き換えても，別名の配列の中身を書き換えても，どちらでも同様に書き換えられるし，また，同様に中身を読み取ることもできる．

あとは実際の例をみていこう．

1# 普通のベクトル．Vector は1次元配列の別の言い方．
2# 添字の範囲は普通に 1:N+4 になるね．
3v = Vector( -2.0:N+1 ) 
4
5# v の添字の範囲を -3 ずらして，o_v という別の名前をつける
6# よって添字の範囲は -2:N+1 になるはずだ．
7o_v = OffsetArray( v, -3 )

204-element OffsetArray(::Vector{Float64}, -2:201) 
  with eltype Float64 with indices -2:201: 
  -2.0
  -1.0
   0.0
  …略…

上の出力を見るとわかるように，配列 o_v の添字(indices)が -2:201 となっていて，通常の 1始まりでないことが明確だね．

ちなみに，以降間違えないように，OffsetArray として添字がずれている配列には接頭語 o_ をつけているので，それをヒントに読み解こう．

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

1function pBC!(o_u)  # 境界条件に従わせる.
2    o_u[-2] = o_u[N-2]
3    o_u[-1] = o_u[N-1]
4    o_u[N]  = o_u[0]
5    o_u[N+1] = o_u[1]
6    return o_u
7end

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

 1function diff1(o_u)
 2    o_v = copy(o_u)
 3    
 4    for k in 0:N-1
 5        o_v[k] = ( o_u[k+1] - o_u[k-1] )/(2Δx)
 6    end
 7    
 8    return pBC!(o_v)
 9end
10
11function diff2(o_u)
12    o_v = copy(o_u)
13    
14    for k in 0:N-1
15        o_v[k] = ( o_u[k+1] -2o_u[k] + o_u[k-1] ) /(Δx^2)
16    end
17    
18    return pBC!(o_v)
19end
20
21diff3(o_u) = diff1(diff2(o_u))

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

1f(o_u) = - o_u .* diff1(o_u) - ϵ^2 * diff3(o_u)

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

1function RK(o_u)
2    f1 = f(o_u)
3    f2 = f(o_u + Δt/2 * f1)
4    f3 = f(o_u + Δt/2 * f2)
5    f4 = f(o_u + Δt * f3)
6    return o_u + Δt * (f1 + 2*f2 + 2*f3 + f4)/6
7end

Zabusky and Kruskal で使われている初期値を設定しよう．

1o_u0 = similar(o_v) # 上で o_v を作っているのでその型を使って．
2
3for k in 0:N-1
4    o_u0[k] = cos(π*k*Δx) # cos の中身は、円周率パイ * k * デルタx 
5end
6
7pBC!(o_u0)

試しに、初期値をプロットしてみる．

1using Plots
2plot( o_u0, marker = :circ )

svg

ふむ、こんな感じか．

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

 1o_u = copy( o_u0 )
 2sq_u = copy( o_u0[0:N-1] )  
 3# 内部点のデータが本質的なので，そこだけ記録する．
 4
 5for i in 1:1000
 6    o_u = RK( o_u )
 7    if i % 200 == 0  
 8    # そのままだとデータが多すぎるので、
 9    # 200ステップおきにデータを格納する
10        sq_u = hcat( sq_u, o_u[0:N-1] )
11    end
12end

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

1sq_u

200×6 Matrix{Float64}:
 1.0       0.859704  0.67817    0.602996    -0.538929   -0.667305
 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.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

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

1plot(sq_u[:,1], marker = :circ )

svg

最後の結果を描いてみる

1plot( sq_u[:,end], marker = :circ )

svg

一応、全部重ねて描いてみる(見にくいのでマーカーは描いてない).

1plot( sq_u )

svg

さらに，これを動画にしてみたのが以下のものだ(→レポート問題)． gif

何が起きているか？

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

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

wikipedia の右下に、同じ計算を motion gif にしたものが掲載されているので見てみよう．

計算結果のチェック: KdV 方程式の保存量をチェックしよう

実は，KdV 方程式には、無限個の保存量があるという，大変変わった特徴があることがわかっている．単純な方から挙げてみると、

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

という感じだ． それぞれ、プログラムしてみよう(一応、入力が普通の配列でも OffsetArray でも大丈夫なように書いてみた).

 1I1(u) = sum(u[1:N]) * Δx
 2# 台形則だ．周期的境界条件なのでこれで正しい
 3
 4I2(u) = sum( u[1:N].^2 ) * Δx
 5# これも台形則だ．周期的境界条件なので以下同文．
 6
 7function I3(u)
 8
 9    r1 = sum( (u.^3)/3 )
10    r2 = sum( (u[k+1] - u[k])^2 for k in 1:N-1 )/(Δx^2)
11    r2 += (u[1]-u[N])^2/(Δx^2) # 面倒だけど、こいつもいれておかないとね…
12    
13    return (r1 - ϵ^2 * r2)*Δx
14end

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

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

6-element Vector{Float64}:
 -3.552713678800501e-17
 -2.1094237467877975e-17
 -1.176836406102666e-16
 -2.042810365310288e-16
 -2.6423307986078724e-16
 -2.486899575160351e-16

ふむ，どうやら保存量 $I_1$ の値はこの初期値だと 0 のようだな．

では $I_2$ はどうか．

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

6-element Vector{Float64}:
 1.0
 1.0001464497535437
 1.0057580637364603
 1.0233388432308645
 1.0279143024252264
 1.0258796813369369

ふむ，$I_2$ の値は本来は 1.0 のようだ．

さて $I_3$ をみてみよう．

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

6-element Vector{Float64}:
 -0.004776495659718536
 -0.004778626109137531
 -0.003959257159688648
 -0.0008801863505688345
 -0.0008322497315909061
 -0.001360199737689527

ふ～む，$I_3$ の値は 0.0048 ぐらいか？ しかもこの計算結果だとあまり保存されてない．

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

計算結果のチェック: 複数のソリトンの挙動を追いかけてみよう

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

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

(疑似) 2-ソリトン解

(数学による大変驚異的な成果であるのだが)高さ $H$, 中心 $x_0$ の 1-ソリトン解は次の式で書ける．

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

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

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

注: 2-ソリトン解は 1-ソリトン解を2つ並べたものとは数式も形状も異なるが、ソリトン同士がある程度離れている場合は 1-ソリトン解 2つでかなりよく近似できる． なにしろ 2-ソリトン解の数式は面倒なので使いたくない. その式を見てみたい人は，wikipedia の該当する部分を見てみよう．

1for k in 0:N-1
2    o_u0[k] = OneSoliton(1, 0.5, k*Δx)+OneSoliton(0.5, 1.2, k*Δx)
3end
4pBC!(o_u0)

204-element OffsetArray(::Vector{Float64}, -2:201) with eltype Float64 with indices -2:201:
 1.034997625577477e-6
 8.597031730922938e-7
 8.006527803097905e-6
 1.0409137026988382e-5
 1.3532723405147744e-5
 ⋮
 1.2460347712918558e-6
 1.034997625577477e-6
 8.597031730922938e-7
 8.006527803097905e-6
 1.0409137026988382e-5

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

1plot( o_u0, marker = :circ )

svg

背の高い波と低い波があることがわかるだろう． そして，左の背の高い波がより速く移動するので，いずれ低い方を追い越すはずだ．その様子が計算結果で見えると良いな，というところだ．

さて，早速計算してみよう

1o_u = copy( o_u0 )
2sq_u = copy( o_u0[0:N-1] )
3
4for i in 1:5000
5    o_u = RK(o_u)
6    if i % 100 == 0
7        sq_u = hcat( sq_u, o_u[0:N-1] )
8    end
9end

ではプロットしてみよう

1plot(sq_u, legend = :no)

svg

これを動画にしてみたものが次のものだ(→レポート問題)． gif

こういう「軌跡」を見るための定番の手法も使ってみよう. 上から見て みる，のだ．

1X = 0:Δx:2-Δx
2T = 0:100Δt:5000Δt
3plot(X, T, sq_u', xlabel = "space x", ylabel = "time t")

svg

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

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

この計算では保存量は？

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

51-element Vector{Float64}:
 0.26019771077247456
 0.26019771077247456
 0.2601977107724746
 0.2601977107724745
 0.2601977107724746
 ⋮
 0.26019771077247444
 0.2601977107724745
 0.26019771077247444
 0.26019771077247444
 0.2601977107724745

ふむ，$I_1$ はやはりよく保存されているな．

では次だ．

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

51-element Vector{Float64}:
 0.1375433423024314
 0.13754627833964014
 0.13754746146211605
 0.1375478920264056
 0.1375475749951093
 ⋮
 0.1375199717151598
 0.13752654005383488
 0.1375316113604702
 0.1375372483263182
 0.1375389858321474

ふ～む，$I_2$ の保存性もなかなか悪くない．

では $I_3$ はどうだろう．

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

51-element Vector{Float64}:
 0.02394673717832975
 0.023947303333121436
 0.02394773044682373
 0.02394787829394399
 0.023947824030209928
 ⋮
 0.023934525788413052
 0.02393779109050128
 0.023940302733468456
 0.023942507565480175
 0.023943729883240412

数字をみるだけなら，これも悪くなさそうだ．

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

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

svg

┌ Warning: No strict ticks found
└ @ PlotUtils c:\julia\PKG\v1.8\packages\PlotUtils\igbcf\src\ticks.jl:191
┌ Warning: No strict ticks found
└ @ PlotUtils c:\julia\PKG\v1.8\packages\PlotUtils\igbcf\src\ticks.jl:191

上下の動きが小さすぎるようで，プロットの警告が出ているし，軸に数字も表示されていない． よく保存しているということなので，まあ嬉しい悲鳴ということだが，念のため数字そのものも見ておこう．

1sq_I1 .- sq_I1[1]

51-element Vector{Float64}:
  0.0
  0.0
  5.551115123125783e-17
 -5.551115123125783e-17
  5.551115123125783e-17
 …略…
 -1.1102230246251565e-16
 -5.551115123125783e-17
 -1.1102230246251565e-16
 -1.1102230246251565e-16
 -5.551115123125783e-17

…まあ見てわかるように，倍精度の数字としてはブレはほぼゼロと言っていいんだろう．

さて，$I_2$, $I_3$ については次のような感じだ．

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

svg

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

svg

そうだなあ，やはり，保存量 $I_1$ は数値的によく保存されているが、続く2つの $I_2$, $I_3$ の保存性はあまり良くない． とくに $I_3$ はもとの値と比べるとパーセントオーダーで誤差が発生しているね．

そして、その誤差の発生タイミングを見ると、どうやら、二つのソリトンが衝突するあたりで結構な数値誤差が発生しているようだということがわかる．

レポート

下記要領でレポートを出してみよう．

  • e-mail にて，
  • 題名を 2024-numerical-analysis-report-09 として,
  • 教官宛(アドレスは web の "TOP" を見よう)に，
  • 自分の学籍番号と名前を必ず書き込んで，
  • 内容はテキストでも良いし，pdf などの電子ファイルを添付しても良いので，

下記の内容を実行して，結果や解析，感想等をレポートとして提出しよう．

  注意
　　近年はセキュリティ上の懸念から，実行形式のプログラムなどをメールに添付するとそのメールそのものの受信を受信側サーバが拒絶したりする． そういうことを避けるため，レポートをメールで提出するときは添付ファイルにそういった懸念のあるファイルが無いようにしよう．

まあ要するに，レポートは pdf ファイルにして，それをメールに添付して送るのが良い ということだと思っておこう．

  1. KdV 方程式の近似解の時間発展の様子(できれば最初の例と，2ソリトンの例の両方)をどうにかして自分で動画にしてみよう．
    なお，そのときに動画を作った方法やプログラムも一緒に提出すること．

  2. KdV 方程式の数値計算を、今使った Runge-Kutta スキームの部分をなにか他の解法に替えてやってみよう．ただし，Euler法でやってみるとひどい目に遭うので，Euler法を試すなら覚悟しておこう．

  1. もともとの境界条件からみて $u_N(t)$ は $u_0(t)$ と同一なので，$u_N(t)$ を独立な存在ではない…ということから，計算対象に含めずに計算対象は $u_{N-1}(t)$ までとするこの形が一番純朴で素直かな． ． ↩︎