A4. 差分法 for 偏微分方程式

Daisuke Furihata (降籏 大介)

grid Photo by Willian Justen de Vasconcellos on Unsplash

偏微分方程式の数値解法: もう少し踏み込む手法の一つ．有限差分法(finite difference method).

これまで何度か登場したが、「微分」を離散近似するシンプルな手法に、差分法と呼ばれるものがある． これは微分を、引き算と割り算で近似するもので、 method of line の回で右辺を離散近似する際でも使われていた．

この差分法を全面的に使って、微分方程式の微分をすべて差分で近似してしまって、偏微分方程式を完全に離散的な式(スキーム)で近似してしまおう、というのが今回紹介する方法だ．

この手法は単純に見えるが、単純なだけに応用は広いし、また、近似の仕方のバリエーションが豊富なため実は奥も深い． そしてなにより、うまくはまれば大変強力な数値解法になる．

というわけで偏微分方程式の数値解法としては主流の方法のひとつなのだが，授業としては既に method of line を説明したことからこの手法はあまりそこから変化したように見えないと思う．そこでこの内容については Appendix の一つとした．

では，この手法についてまずは入門的な内容で様子をみてみよう．

対象は KdV 方程式で

method of line の回で紹介した KdV 方程式として, パラメータ $\epsilon$ をもつ次の式

\[ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \epsilon^2 \frac{\partial^3 u}{\partial x^3} = 0 \]

を考えることにしよう．

差分法を使ってみる

差分法は微分を差分、つまり、Taylor 展開で近似するものである． 近似できていれば良いのだ！ と考えれば、こんなに自由な方法もないわけで、まあやりたい放題といって良い． 自由であるということと、良い近似解法ができるかどうかは全く別だが．

そこで、たとえば今回は、あちこち平均したり差分したりして作った以下の離散式(数値スキーム)、をもとの KdV 方程式の近似式として考えよう． ただし、$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} + \left( \frac{u^{+} + u}{2} \right) \delta_k^{(1)} \left( \frac{u^{+} + u}{2} \right) + \epsilon^2 \delta_k^{(3)} \left( \frac{u^{+} + u}{2} \right) = 0 \label{eq:scheme} \]

ただし、$n = 0,1,2, \cdots$, $k = 0, 1, \cdots, N-1$ で、離散境界条件は以前と同じように周期的境界条件を離散化した

\[ \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. \label{eq:scheme-bc} \]

としよう．

早速計算してみよう！

では早速計算へ向けて準備を始めよう．

まずは method of line の回で説明したのと同様に， wikipedia でも引用されている Zabusky and Kruskal の1965年の古い論文 で使われているパラメータ等にあわせよう．

1
2
3
4
5

const ϵ,L = 0.022, 2

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

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

1

using OffsetArrays

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

1
2
3
4
5
6
7

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

初期値等も method of line の回で説明したのと同様に定義しよう． まず，初期値を作ってしまう．

1
2
3
4
5

# 通常の array. 添字は 1 から．Vector は１次元配列の別名.
u0 = Vector( -2.0:N+1 ) 

# u0 の添字を -3 ずらす．つまり，-2 からになるね．
o_u0 = OffsetArray( u0, -3 )

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

1
2
3
4

# 境界内部は定義通り計算して代入．
o_u0[0:N-1] = [ cos(π*k*Δx) for k in 0:N-1 ]  

pBC!(o_u0)   # 境界条件を満たすように, 外部仮想点の値を設定．

204-element OffsetArray(::Vector{Float64}, -2:201) with eltype Float64 with indices -2:201:
 0.9980267284282716
 0.9995065603657316
 1.0
 0.9995065603657316
 0.9980267284282716
 0.99556196460308
 0.9921147013144779
 ⋮
 0.9921147013144778
 0.99556196460308
 0.9980267284282716
 0.9995065603657316
 1.0
 0.9995065603657316

これまた method of line の回で説明したのと同様，一階から三階までの差分作用素を関数として定義しておこう．

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

function diff1(o_u)
    o_v = similar(o_u)
    o_v[0:N-1] = [( o_u[k+1] - o_u[k-1] )/(2Δx) for k in 0:N-1 ] 
    return pBC!(o_v)
end

function diff2(o_u)
    o_v = similar(o_u)
    o_v[0:N-1] = [ ( o_u[k+1] -2o_u[k] + o_u[k-1] ) /(Δx^2) for k in 0:N-1 ]  
    return pBC!(o_v)
end

diff3(o_u) = diff1(diff2(o_u))

数値スキームそのものを関数にしてしまおう

さて，スキームの式 $(\ref{eq:scheme})$ は見てわかるように，新しい値 $\{u_k^{+} \}_{k=0}^{N-1}$ に対して非線形(連立)方程式の形をしている．

そのため，これを解くために NLsolve package を使うとしよう．

それには $0 = 式$ という形に式を直しておかないといけない．まあ今回は既にそうなっているが． あと，式の形を見るとわかるだろうが，ついでに式全体に $\Delta t$ をかけておくと，式全体のオーダーが $u$ と同じオーダーになるのでそうしておこう．

1
2
3
4

function f(o_up, o_u)
    o_av = ( o_up + o_u )/2
    return o_up - o_u + Δt*( o_av .* diff1(o_av) + ϵ^2 * diff3(o_av) )
end

そして、この式がゼロになるような解 up を「探す」ことで時間ステップが進んだ値である up を求めるわけだ(その計算自体は NLsolve package がやってくれるね)．

さて，あとは具体的にどの変数を NLsolve に渡して解いてもらうかを考えよう．

それにはおおざっぱに２つのやり方がある．

  1. 最初の方法は，「すべてのデータを渡し，すべての条件を解いてすべての解を求めてもらう」原理的なものだ．

    つまり，外部仮想点も含めたデータ $\{u_k \}_{k=-2}^{N+1}$ を渡し， スキームの式 $(\ref{eq:scheme})$ と 境界条件式 $(\ref{eq:scheme-bc})$ の両方の式も渡してこれを満たすように 解 $\{ u_k^{+} \}_{k=-2}^{N+1}$ を計算してもらう，という手だ．

    これは一見悪くない手なのだが，実は良くない． というのも， スキームの式 $(\ref{eq:scheme})$ は微分方程式(の近似)で， 境界条件式 $(\ref{eq:scheme-bc})$ は代数方程式であり，数学的にだいぶ性質が異なるのだ． これらを「数値的に同時に解こうとする」のは「DAE 問題を解かないといけない」という，実はちゃんとした難しい問題1で，NLsolve package が答えを返さないという可能性も結構高い．

  2. ２つ目の方法は，境界条件式 $(\ref{eq:scheme-bc})$ を先にスキームの式に反映させてしまい，そのスキーム部分だけを解いてもらう，という方法だ．

    つまり，外部仮想点も含めないデータ $\{u_k \}_{k=-0}^{N-1}$ を渡し， 境界条件を考慮した形でのスキームの式 を渡してこれを満たすように 解 $\{ u_k^{+} \}_{k=0}^{N-1}$ を計算してもらう，という手だ．

    これはいかにも恣意的な解法のように見えるが，1. が抱える「DAE を解かないといけない」という問題を回避しており，実はこちらが推奨される方法となる．

というわけで，上に解説したように，今回は境界条件を先にスキームに反映させ，その変形した形のスキームを解いてもらうことにしよう． 幸い，差分作用素の関数は境界条件を反映するように作ったので，現状でスキームは「そうなっている」．

なので，あとは 外部仮想点も含めないデータ $\{u_k \}_{k=-0}^{N-1}$ を渡し， 解 $\{ u_k^{+} \}_{k=0}^{N-1}$ を解いてもらい，結果を 新値 $\{ u_k^{+} \}_{k=-2}^{N+1}$ に拡張する，ということを繰り返せばよい．

そのために，内部点のみのデータと，外部仮想点も含むデータを相互に変換できるように関数を作っておこう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# 内部点を取り出す．
inside(u) = u[0:N-1] 

# 内部点データから，外部仮想点も含めたデータへ拡張．
function expand(u)
    o_u = similar(o_u0) # 初期値の offsetarray o_u0 と同じ型．
    o_u[0:N-1] = u #  はみ出ていない中身は本質的に同じ．添字が一つずれてるけどね．
    pBC!(o_u)
    return o_u
end

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

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

function f_solve(up, u)  
# NLsolve の為、入力値の up も u も標準の配列でないといけない．

    o_up, o_u = expand(up), expand(u)
    # スキームは Offset array 用なので変数を拡張．

    o_f = f(o_up,o_u)
    # スキームの「残差(ベクトル)」を計算．

    return inside( o_f ) 
    # NLsolve の為、出力も標準の配列でないといけない．
end

これで良い．

あとはこの f_solve 関数がゼロになるように数値解 up を探せ！ と NLvole に命令すれば良いことになる． まずいつものように 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

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

1

u1 = nls(f_solve, inside(o_u0), ini= inside(o_u0) )

200-element Vector{Float64}:
0.9999951869757729
0.9995999059182815
0.9982178772929537
0.9958501739167067
0.9924988433939419
⋮
0.9872009729451615
0.9917214646657009
0.9952644297153576
0.997826087542708
0.9994036223443052

  環境によってはこの計算時に StarkOverflowError というエラーが出て実行できないことがある． これは「関数の多重呼び出しの深さが深すぎた」ことにより起きるエラーで，まあ簡単に言うと，あまり考えずに関数を重ねていくようなプログラムを組むとこうなるとでも言えば良いかな．

解決するにはなるべく関数呼び出しがあまり多重にならないようなプログラムに書き直すのが王道で，今回は inside や expand を使わずに済むように 書き換えるというあたりになるだろうか． ただまあそれも面倒なので，今回の場合は NLsolve で使うアルゴリズムをよりシンプルなものに変えてもらうようにオプションで指定するのが良い． 具体的には，method = :anderson というオプションを使うように nls 関数を定義しなおせば良いだろう．次のようにすれば良いだろう．

1
2
3
4
5
6
7

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, method = :anderson).zero
    end
end

さて，話を戻し，計算結果を図にしてみて、変でないことをたしかめておこう．

1
2

using Plots
plot( u1 )

ふむ、こんな感じか． 破綻等はしていないようで、まずは一安心．

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

その前に, 進度をリアルタイムで見られるように

このスキームは非線形方程式を解くので、時間がかかることが想定される． そこで、計算がどれくらい進んでいるのか、リアルタイムで見えるようにしよう．それには、ProgressMeter というパッケージを使えば良い.

インストールしてない人は，いつものように

1
2

using Pkg
Pkg.add("PrograessMeter")

としてインストールしておこう．

これを使うと、@showpprogress 時間間隔(秒) メッセージ などを for 文の頭につけるだけで、その for 文の実行状況と、残り推定時間を指定した時間間隔毎に更新して教えてくれる． 使い方はこれから見る例で十分だろう．

ではまず，使用宣言をしよう．

1

using ProgressMeter

これでよし．では計算しよう．

ただし、次の計算は関数 nls に method = :anderson オプションを付けていない場合は計算時間がけっこうかかる． そうした場合は N や for 文の数を小さくするなど、現実的な計算時間になるようにパラメータを調整するか，このオプションを付け直してやり直そう.

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

u = inside(o_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  
   # 2step 目から、推測値を「過去と今から計算」して入力．
   # こうすると計算が少し速い．

   if i % 100 == 0  # 100ステップおきにデータを格納するとして
        sq_u = hcat(sq_u, u)
    end
end

Progress: 100%|████████| Time: 0:00:00 ← method = :anderson オプション有りの場合

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

1

sq_u

200×11 Array{Float64,2}:
1.0       0.956153  0.859704  0.757938  …  -0.538929   -0.66376    -0.667301
0.999507  0.964282  0.871752  0.770684     -0.576519   -0.630957   -0.602904
0.998027  0.971663  0.88342   0.783294     -0.559205   -0.553665   -0.497812
0.995562  0.978272  0.89469   0.79573      -0.479927   -0.417903   -0.336871
0.992115  0.984084  0.905545  0.807932     -0.332199   -0.215148   -0.110155
⋮                                       ⋱                           ⋮       
0.987688  0.905118  0.794297  0.692256  …   0.551663   -0.133731   -0.487255
0.992115  0.916636  0.80802   0.705686      0.220001   -0.349151   -0.59546 
0.995562  0.92752   0.821436  0.718959     -0.0644842  -0.500781   -0.660936
0.998027  0.937749  0.834531  0.732083     -0.287025   -0.600866   -0.69433 
0.999507  0.947301  0.847292  0.745072     -0.442895   -0.652234   -0.695196

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

1
2

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


1
2

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


1
2

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

以前同様、ソリトンへの分解がきちんと見えている．

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

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

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

method of line の回でも書いたように、KdV 方程式には、無限個の保存量があることがわかっていて、単純な方から挙げてみると、

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

という感じだ． それぞれ method of line の回と同様にプログラムしてみよう．

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
2

len = size(sq_u)[2]
sq_I1 = [ I1(sq_u[:,k]) for k in 1:len ]

11-element Vector{Float64}:
 -3.552713678800501e-17
  4.107825191113079e-17
 -1.3322676295501878e-17
 -1.9206858326015208e-16
 -1.9095836023552692e-16
  3.6415315207705137e-16
  1.821875983409882e-15
  1.8807178037150153e-15
  1.9345636204093354e-15
  1.9118040484045194e-15
  1.9328982858723976e-15

ふむ，$I_1$ は値 0 としてほぼ保たれているようだ．

では次に $I_2$ はどうかな．

1

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

11-element Vector{Float64}:
 1.0
 1.0000265358127824
 1.00014645856071
 1.000823350458995
 1.0057579344512129
 1.0157310084577404
 1.0233389542936338
 1.026994465186611
 1.0279146284586456
 1.0272991063204429
 1.0258802772889797

ふ～む，$I_2$ には 3% 弱のズレが見られるか．

では $I_3$ はどうかな．

1

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

11-element Vector{Float64}:
 -0.004776495659718536
 -0.004776884780748252
 -0.00477862419108704
 -0.004756513550836625
 -0.003959081321684294
 -0.0020069956259172626
 -0.0008799771146812674
 -0.0006626322917696115
 -0.0008319409212678864
 -0.0011006884030202002
 -0.001359871995866726

う～ん… $I_3$ はそれなりにずれているようだ．

とりあえず、この数値解法だとこの 3つの保存量はおおよそこのような感じになることがわかる．

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

先にも書いたように、ソリトンはピーク値が大きい(背が高い)ほうが速く移動する．

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

(疑似) 2-ソリトン解

先にも示したように、高さ $H$, 中心 $x_0$ の 1-ソリトン解は次の式で書ける．

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

そこでこれまた method of line の回と同様、少し離してこれを二つ配置することで 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

o_u0[0:N-1] = [ OneSoliton(1, 0.5, k*Δx) + OneSoliton(0.5, 1.2, k*Δx) for k in 0:N-1]
pBC!(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

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

1

plot( o_u0 )

問題なさそうだ．

では早速計算してみよう． これまた、計算時間が膨大になりそうなときは、現実的な計算時間になるようにパラメータを調整してやり直そう.

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

Δt = 0.001

u = inside( o_u0 )
sq_u = copy( u )

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

@showprogress for i in 2:5000  

   u, u_old = nls( f_solve, u, ini = 2*u - u_old ), u
   # 2step 目から、推測値を「過去と今から計算」して入力．
    
   if i % 50 == 0  # 50ステップおきにデータを格納するとして
        sq_u = hcat(sq_u, u)
    end
end

Progress: 100%|██████████| Time: 0:00:02 ← method = :anderson の場合

ではプロットしてみよう

1

plot( sq_u, legend = :no )

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

1
2
3

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

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

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

この計算では保存量は？

計算してみよう．

1
2

len = size(sq_u)[2]
sq_I1 = [ I1(sq_u[:,k]) for k in 1:len ]

101-element Vector{Float64}:
 0.26019771077247456
 0.2601977107724745
 0.2601977107724746
 0.2601977107724746
 0.26019771077247467
 ⋮
 0.26019771077247433
 0.2601977107724744
 0.26019771077247433
 0.2601977107724744
 0.26019771077247444

ふむ，良いね．$I_1$ の保存性は高いようだ．

では $I_2$ はどうか．

1

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

101-element Vector{Float64}:
 0.1375433423024314
 0.13754489027309144
 0.13754628165489616
 0.13754704620105346
 0.13754746630880768
 ⋮
 0.13753164070421137
 0.137534597309656
 0.13753726243951353
 0.13753865541386553
 0.1375390191334265

これも一見では悪くないな．では $I_3$ は？

1

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

101-element Vector{Float64}:
 0.02394673717832975
 0.023946911504337493
 0.02394730410837376
 0.02394757752431669
 0.02394773168947462
 ⋮
 0.02394031115024523
 0.02394147091124976
 0.023942522933476144
 0.023943285323210457
 0.023943739417863732

これもまあ悪くなさそうに見える．

では、それぞれグラフにしてみよう． method of line の回同様、変化が小さいので、初期の値を引いてそこからのズレをみると良いだろう．

1

plot(sq_I1 .- sq_I1[1] , marker = :auto, legend = :no ) 


┌ Warning: No strict ticks found
└ @ PlotUtils c:\julia\PKG\v1.6\packages\PlotUtils\VgXdq\src\ticks.jl:294
┌ Warning: No strict ticks found
└ @ PlotUtils c:\julia\PKG\v1.6\packages\PlotUtils\VgXdq\src\ticks.jl:294

1

plot( sq_I2 .- sq_I2[1], marker = :auto, legend = :no )


1

plot( sq_I3 .- sq_I3[1], marker = :auto, legend = :no )

保存量 $I_1$ は数値的に大変よく保存されているね．ほぼ完璧だ．

他の 2つの保存量 $I_2$, $I_3$ の保存性はソリトンの衝突のあたりでやっぱりあまり良くない． このあたりは method of line の回の method of line + Runge-Kutta 法による結果とよく似ているね．

まあ，結果としては差分法も悪くないね．

非線形連立方程式の形の差分スキームだと少し計算時間がかかるけど，本当に計算時間が気になるなら予測子修正子法を適用するという手もあるしね．

  1. こうした微分方程式と代数方程式の連立したものを微分代数方程式 DAE (Differential-Algebraic Equations) と呼び，数値的に解くのが難しいことが知られている．それ自体が研究対象なのだ．厄介なのだ． ↩︎