grid Photo by Willian Justen de Vasconcellos on Unsplash

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

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

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

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

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

対象は KdV 方程式で

前回紹介した 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} \]

としよう．

早速計算してみよう！

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

まずは前回同様， 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

初期値等も前回同様に定義しよう．まず，初期値を作ってしまう．

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(::Array{Float64,1}, -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

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

 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 に渡して解いてもらうかを考えよう．

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

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

つまり，外部仮想点も含めたデータ $\{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 問題を解かないといけない」という，実はちゃんとした難しい問題¹で，NLsolve package が答えを返さないという可能性も結構高い．
２つ目の方法は，境界条件式 $(\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 Array{Float64,1}:
0.9999951869757729
0.9995999059182815
0.9982178772929537
0.9958501739167067
0.9924988433939419
⋮
0.9872009729451615
0.9917214646657009
0.9952644297153576
0.997826087542708
0.9994036223443052

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

1
2


using Plots
plot( u1 )

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

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

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

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

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

1
2


using Pkg
Pkg.add("PrograessMeter")

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

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

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

1

using ProgressMeter

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

ただし、次の計算は遅いマシンの場合、計算時間がたいへんかかる (速めの PC だと 1分程度なんだけどね)．そうした場合は N や for 文の数を小さくするなど、現実的な計算時間になるようにパラメータを調整してやり直そう.

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


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
   # 普通の PC だと 3分ぐらいかかるか？ 

   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:58

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

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 )

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

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

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

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

前回も書いたように、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. \]

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

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 Array{Float64,1}:
-3.552713678800501e-17
-1.6542323066914832e-16
-6.494804694057166e-16
-1.5165646516379638e-15
-7.760458942129844e-16
-7.327471962526034e-16
-8.459899447643693e-16
4.894695759816159e-16
2.473576898864849e-15
6.257216966787382e-15
9.571232695293475e-15

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

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

1

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

11-element Array{Float64,1}:
1.0
1.000026531922912
1.0001464491334295
1.0008233353079634
1.0057579073883398
1.0157309590265462
1.023338890469956
1.0269943893725628
1.0279145401745813
1.0272990060416662
1.025880166582934

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

では $I_3$ はどうかな．

1

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

11-element Array{Float64,1}:
-0.004776495659718536
-0.004776884877922688
-0.004778624463676382
-0.004756513947220731
-0.003959083827698589
-0.0020070016825051696
-0.0008799844342889429
-0.0006626407905056908
-0.0008319507341631294
-0.00110069924584856
-0.0013598837884592286

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

とりあえず、この数値解法だとこの 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

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(::Array{Float64,1}, -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
18
19


Δt = 0.002 # 計算時間がかかるので、時間ステップ間隔を大きくしてしまえ．

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

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

@showprogress for i in 2:2500  
  # 遅いマシンだと 5分以上かかるかも ．
  # 遅くて辛い場合は適当にパラメータを変えよう．

   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:02:25

ではプロットしてみよう

1

plot( sq_u, legend = :no )

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

1
2
3


X = 0:Δx:2-Δx         # 真面目に x の範囲を書き…
T = 0:50Δt:2500Δ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 ]

51-element Array{Float64,1}:
 0.26019771077247456
 0.2601977107724769
 0.26019771077247916
 0.26019771077248144
 0.26019771077248377
 ⋮
 0.2601977107725766
 0.260197710772579
 0.26019771077258164
 0.26019771077258413
 0.26019771077258635

ふむ，悪くない．$I_1$ の保存性は高いようだ．

では $I_2$ はどうか．

1

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

51-element Array{Float64,1}:
 0.1375433423024314
 0.1375462801735097
 0.13754746289081918
 0.13754789130759676
 0.13754757792031222
 ⋮
 0.1375199490514289
 0.13752652170133822
 0.13753160374419313
 0.13753718483703747
 0.1375389787618209

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

1

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

51-element Array{Float64,1}:
 0.023946737178329753
 0.023947302658895858
 0.02394772972563777
 0.023947878264761756
 0.023947821775390184
 ⋮
 0.023934507188682976
 0.023937780539564697
 0.02394029388482077
 0.023942517037247806
 0.023943730562586713

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

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

1

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

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 + Runge-Kutta 法による結果とよく似ているね．

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

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

レポート

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

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

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

4つ目の保存量チェック.
wikipedia の「保存量」項を参照して4つ目の保存量 $I_4$ を調べて、この保存量 $I_4$ (の離散近似値)が数値計算でどう変化するか、前回(method of line)や今回(差分法)の計算方法に対して調べてみよ．グラフも描こう．
KdV 方程式の数値計算に際し，今回の授業では差分スキーム $(\ref{eq:scheme})$ を NLsolve package で直接解いたがこれは時間がかかる．

そこで，自分で適当に Euler スキームを作ってそれを予測子とし，差分スキーム $(\ref{eq:scheme})$ を適当に移項して修正子とする形で予測子修正子法で今回の差分法を高速に解けないか試みてみよ．

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

数値計算法基礎 (2020)

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