偏微分方程式: 差分法

偏微分方程式の数値解法: もう少し踏み込む手法の一つ.有限差分法(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 \]

ただし、$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. \]

としよう.

早速計算してみよう!

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

まずは 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
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!(u)  # 境界条件に従わせる.
    u[-1] = u[N-1]
    u[-2] = u[N-2]
    u[N]  = u[0]
    u[N+1] = u[1]
    return u
end

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

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)   # 境界条件を満たすようにする

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.9876883405951378 0.9822872507286887 0.9759167619387474 0.9685831611286311 0.9602936856769431 0.9510565162951535 ⋮ 0.9510565162951535 0.960293685676943 0.9685831611286312 0.9759167619387474 0.9822872507286887 0.9876883405951377 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(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))

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

(すぐ後で書くように NLpackage で使う形式にするために) $0 = 式$ となるように移項して、この式部分を記述する格好にする. ただし、全体に $\Delta t$ をかけておくと良いかなあ. 式全体のオーダーを $u$ に合わせると、なにかと見易いのだ.

1
2
3
4
function f(up, u)
    av = (up+u)/2
    return up - u + Δt * ( av .* diff1(av) + ϵ^2 * diff3(av) )
end

そして、この式がゼロになるような解 up を「探す」ことで時間ステップが進んだ値である up を求めるわけだが、それには以前も使った NLsolve package が使える(お手軽だね!).

しかし、NLsolve package は通常の配列しか知らないので、NLsolve package へそうした形で情報を渡せるように、通常の、外へはみ出てない配列と offsetarray を相互に変換する関数を NLsolve 用に作っておこう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# Offset array から標準の array への変換.こっちは簡単.
toArray(u) = u[0:N-1] 

# 標準の array から Offset array への変換.境界データが **pBC!** を使って付与される.
function toOffset(u)
    ou = similar(u0) # 初期値 u0 と同じ型だしね.
    ou[0:N-1] = u #  はみ出ていない中身は本質的に同じ.添字が一つずれてるけどね.
    pBC!(ou)
    return ou
end 

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
function f_solve(up_in, u_in)  
# NLsolve の為、入力値の up_in も u_in も標準の配列でないといけない.

    up, u = toOffset(up_in), toOffset(u_in)
    # スキームは Offset array を必要とするのでそれを作っておいて、

    out_f = f(up,u)
    # スキームの「値(ベクトル)」を計算する.これがゼロベクトルになって欲しい…

    return toArray( out_f ) 

# NLsolve の為、出力も標準の配列でないといけない.
end

これで良い.

あとはこの f_solve 関数がゼロになるように数値解 up_in を探せ! と 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, toArray(u0), ini= toArray(u0) )

200-element Array{Float64,1}: 0.9999951869757729 0.9995999059182815 0.9982178772929537 0.9958501739167067 0.9924988433939419 0.9881669072393687 0.9828583590284564 0.9765781615725981 0.9693322431175374 0.9611274925643408 0.9519717537132879 0.9418738185322493 0.9308434194520916 ⋮ 0.9287042830549835 0.9398818166265805 0.9501347776426452 0.9594528192246948 0.9678265062156535 0.9752473254252416 0.9817076950057261 0.9872009729451615 0.9917214646657009 0.9952644297153576 0.997826087542708 0.9994036223443052

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

1
2
using Plots
plot( u1 )

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

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

その前に

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

今の JuliaBox にはデフォルトでインストールされているので、JuliaBox 利用者はそのまま使える. 自分の PC などでは、 Pkg.add コマンドでインストールしておこう.

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

1
using ProgressMeter

これでよし.では計算しよう. ただし、次の計算は juliabox の場合、サーバの混み具合によっては計算時間がたいへんかかる (通常の PC だと 1分程度なんだけどね). そうした場合は N や for 文の数を小さくするなど、現実的な計算時間になるようにパラメータを調整してやり直そう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
u = toArray(u0)   
# 以降は Offsetarray は関数 f_solve の中などに隠蔽され、見かけ上は登場しない.
# そこで、以降は標準の配列だけが表に出る格好で書こう.

sq_u = copy( u )

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

@showprogress for i in 2:1000
   # 環境にもよるが、3分ぐらいかかるか? 

   global u, u_old, sq_u
   # 関数の中などと同じく、for 文中では変数はデフォルトでは局所変数なので、
   # for 文の外と共通な変数はこうして宣言しておく.

   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:03:01

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

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.989074 0.915966 0.81983 … -0.100764 0.0718194 0.196364 0.982287 0.993214 0.925932 0.831366 0.221209 0.44255 0.57748 0.975917 0.996478 0.935423 0.842516 0.631325 0.881213 1.00676 0.968583 0.998838 0.944416 0.853315 1.08832 1.32207 1.41053 0.960294 1.00027 0.952887 0.86386 1.50778 1.66292 1.68633 0.951057 1.00073 0.960811 0.8743 … 1.76394 1.78613 1.73734 0.940881 1.00021 0.968161 0.884792 1.76168 1.64679 1.54699 0.929776 0.998671 0.974907 0.895442 1.49928 1.29714 1.18418 ⋮ ⋱ ⋮
0.929776 0.808463 0.690631 0.594284 0.412862 1.14712 1.32352 0.940881 0.823853 0.70618 0.608637 0.740933 1.32761 1.14066 0.951057 0.838747 0.721498 0.622884 … 1.04426 1.32333 0.835582 0.960294 0.853127 0.736578 0.63702 1.25063 1.13889 0.489555 0.968583 0.866973 0.751409 0.651037 1.29351 0.829242 0.160083 0.975917 0.880266 0.76598 0.664923 1.15866 0.477859 -0.11554 0.982287 0.892988 0.780279 0.678666 0.886113 0.143983 -0.33151 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
plot(sq_u[:,end]) # 最後の結果を描いてみる

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

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

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

前回も記したが、 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\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

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

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

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

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.759363585357794e-5 2.2873142538428366e-5 2.9736909454011856e-5 3.866032231487881e-5 5.0261403915521214e-5 6.534360397797847e-5 8.495142023661298e-5 0.00011044269472890136 ⋮ 4.567357628844041e-6 3.7938001965000576e-6 3.1512570757272405e-6 2.617539027436925e-6 2.1742148533654365e-6 1.8059749345110413e-6 1.5001025568467903e-6 1.2460347712918558e-6 1.034997625577477e-6 8.597031730922938e-7 8.006527803097905e-6 1.0409137026988382e-5

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

1
plot( toArray(u0) )

では早速計算してみよう. これまた、計算時間が膨大になりそうなときは、現実的な計算時間になるようにパラメータを調整してやり直そう.
(通常の PC だと 3分ぐらいなんだけどね…)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Δt = 0.002 # 計算時間がかかるので、時間ステップ間隔を大きくしてしまえ.

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

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

@showprogress for i in 2:2500  
  # 環境にもよるが、7分30秒ぐらいかかるか? 
  # 授業時は 2:1000 ぐらいにしておいた方が良さ気?    

   global u, u_old, sq_u
   # 関数の中などと同じく、for 文中では変数はデフォルトでは局所変数なので、
   # for 文の外と共通な変数はこうして宣言しておく.

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

ではプロットしてみよう

1
plot(sq_u, legend = :no)

ちょっとわかりにくいので上から見よう.計算パラメータが粗いので、少し図が「ガタついている」のはしかたないね.

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

26-element Array{Float64,1}: 0.26019771077247456 0.26019771077247916 0.26019771077248377 0.2601977107724882 0.2601977107724925 0.2601977107724967 0.2601977107725012 0.2601977107725054 0.26019771077251025 0.26019771077251524 0.26019771077252024 0.26019771077252524 0.2601977107725302 0.2601977107725351 0.26019771077254 0.2601977107725438 0.2601977107725466 0.26019771077254916 0.2601977107725523 0.2601977107725567 0.26019771077256176 0.2601977107725667 0.26019771077257164 0.2601977107725766 0.26019771077258164 0.26019771077258635

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

26-element Array{Float64,1}: 0.1375433423024314 0.13754746289081918 0.13754757792031222 0.13754800010841584 0.1375473997884613 0.1375459116043926 0.13754312760397874 0.13754108512470595 0.13753662493051158 0.13752982152025198 0.13751283427995384 0.13748927360288618 0.1374491452680927 0.13739107833231085 0.13731763407288317 0.13724532552992097 0.13720489908797498 0.1372183805008575 0.1372768951769982 0.1373535773971914 0.13741981885018362 0.13746843550698917 0.13750192038148254 0.1375199490514289 0.13753160374419313 0.1375389787618209

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

26-element Array{Float64,1}: 0.023946737178329753 0.02394772972563777 0.023947821775390184 0.023947879969660447 0.02394767751822008 0.023946974701341283 0.023945959220938437 0.023944807186276323 0.02394249690800102 0.02393864497032242 0.023930684712703005 0.02391872864971971 0.023898387745791085 0.023868112543267622 0.0238293417365997 0.023790750152388875 0.02376878475959214 0.02377581021005925 0.023807560942506575 0.023848211703021552 0.02388302405037638 0.023908169813063953 0.023925042259373985 0.023934507188682976 0.02394029388482077 0.023943730562586713

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

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

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

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

保存量 $I_1$ は数値的によく保存されているが、続く2つの $I_2$, $I_3$ の保存性はあまり良くないことや二つのソリトンが衝突するあたりで数値誤差が発生している様子など、先週の method of line + Runge-Kutta 法による結果とよく似ていることが確認できるだろう.

Report No.10 差分法で計算

  1. 4つ目の保存量チェック
    wikipedia の「保存量」項 を参照して4つ目の保存量 $I_4$ を調べて、これ(の近似値)が数値計算でどう変化するか、前回(method of line)や今回(差分法)の計算方法に対して調べてみよ.グラフも描こう.

  2. 動画チャレンジ
    上の方法で計算した KdV 方程式の数値解を動画にしよう.

  3. 上の差分法スキームを、もっと高速に計算できないだろうか
    いろいろ考えて試してみて、それぞれの計算時間を比較しよう.
    ただし、計算結果があまり崩れるのも意味がなさそうなので、上の3つの保存量チェックぐらいはしておこう.