偏微分方程式: 有限要素法

有限要素法(FEM: Finite Element Method)について

差分法の類は空間の次元が上がっていくと「どう離散化を定義するか」という問題に直面することになる(空間次元が 1次元だと実感しにくいが)． そこで、次元やメッシュの歪みに強い、汎用性の高い方法として有限要素法(finite element method)を扱ってみよう．

まあまずは、FEM についての簡単な手書きメモ(pdf) を読んでみよう．

FEM を使って実際に問題を解いてみよう

KdV 方程式を対象として

ここでもやはり、Koreweg de Vries 方程式として, パラメータ $\epsilon$ をもつ次の式

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

を考えることにしよう．

有限要素法による数値スキームを作ろう

詳細は別資料に記載してあるが、ここでは解を空間方向には区分的一次多項式(つまり折れ線グラフだ)で近似することにして、それぞれのポイントの高さを集めたベクトルを $u$ として、次のような近似式を得る. 時間方向は(まだ)いじってないためこれはいわゆる半離散スキームである．

\[ \Phi \frac{d u}{dt} = - \Psi(u) u - \epsilon^2 D_1 \Phi^{-1} D_2 u \]

ただし、$\Phi, \Psi(u), D_1, D_2$ はそれぞれ行列である．

この半離散数値スキームによってベクトル $u$ に対する時間微分 $du/dt$ が得られるので、 以前紹介した method of line を採用するならば，あとは Runge-Kutta 法などを使えばよいだけだ． 今回はそうしてみよう．

早速計算してみよう！

では早速計算してみよう．常微分方程式の解法には 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
2
3
4
5
6
7

const ϵ,L = 0.022, 2

N = 200
Δx = L/N
Δt = 0.0001  
# やってみるとわかるが、今回のスキームはちょっと繊細で
# Δt を細かくしないとうまく動かない．

Zabusky and Kruskal https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.15.240 で使われている初期値を設定しよう．

1

u0 = [ cos(π*k*Δx) for k in 0:N-1 ]

200-element Array{Float64,1}: 1.0 0.9995065603657316 0.9980267284282716 ⋮ 0.9980267284282716 0.9995065603657316

念のためにプロットして確かめよう(なるべく頻繁にいろいろ確認すると，プログラムのミスに悩む時間が減るぞ)

1
2

using Plots
plot(u0)

ふむ、いつもの初期値だな．ということで確認がとれたので次へ進もう． 次は、スキーム中に出現する行列などを作るとよさそうだ．

注: 下記の SparseArrays package は「すかすか = ゼロ値な要素が多い」行列やベクトルを扱うのに適したパッケージで，Julia では標準装備のものだ．適した問題にはメモリも計算量も少なくて済むことが多いので積極的に使っていこう．

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

using LinearAlgebra
# 線形代数の少し突っ込んだ機能を使うかもしれないのでこの宣言をしておく

using SparseArrays
# すかすかな行列やベクトルを扱うときはメモリや計算が無駄にならないように


v = ones(N-1)
Φ = (2/3)*I + (1/6)* diagm(1 => v, -1 => v)
Φ[1,N] = Φ[N,1] = 1/6
Φ = sparse( Δx * Φ )
# すかすかな行列やベクトルは sparse 扱いをしておく．
# メモリや計算量が少なくて済むぞ．

200×200 SparseMatrixCSC{Float64,Int64} with 600 stored entries: [1 , 1] = 0.00666667 [2 , 1] = 0.00166667 [200, 1] = 0.00166667 [1 , 2] = 0.00166667 [2 , 2] = 0.00666667 [3 , 2] = 0.00166667 [2 , 3] = 0.00166667 [3 , 3] = 0.00666667 [4 , 3] = 0.00166667 [3 , 4] = 0.00166667 [4 , 4] = 0.00666667 [5 , 4] = 0.00166667 ⋮ [196, 197] = 0.00166667 [197, 197] = 0.00666667 [198, 197] = 0.00166667 [197, 198] = 0.00166667 [198, 198] = 0.00666667 [199, 198] = 0.00166667 [198, 199] = 0.00166667 [199, 199] = 0.00666667 [200, 199] = 0.00166667 [1 , 200] = 0.00166667 [199, 200] = 0.00166667 [200, 200] = 0.00666667

注* I は「大きさが自由な」単位行列．このように式の他の項から大きさが自動的に定まるようなケースではこう書くだけで使えて便利だ． 大きさを定めた単位行列を作りたいという場合(昔あった eye という関数は廃止された)は，そうだなあ，例えば

1

Matrix{Float64}(I, 3, 3)

とか

1

diagm(0 => ones(3))

とすると $3 \times 3$ の(実数値の)単位行列が作れるぞ．

さて話を戻して，次に、中心差分を作るとしよう．

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

function cyclic_diff_center(u)
    v = similar(u)
    n = length(u)
    v[1] = u[2] - u[n]
    v[n] = u[1] - u[n-1]
    for k in 2:n-1
        v[k] = u[k+1] - u[k-1]
    end
    return v
end

function Ψ(u)
    vc = (1/3) * cyclic_diff_center(u)
    vf = (1/6) * diff(u)

    corner = (1/6) * ( u[1] - u[end] )

    M = diagm(0 => vc, 1 => vf, -1 => vf)
    M[1,end] = M[end,1] = corner

    return M
end

Ψ (generic function with 1 method)

1
2
3
4
5

v = (1/2) * ones(N-1)
D1 = diagm( -1 => v, 1 => -v )
D1[1,end] = 1/2
D1[end,1] = -1/2
D1 = sparse( D1 )

200×200 SparseMatrixCSC{Float64,Int64} with 400 stored entries: [2 , 1] = 0.5 [200, 1] = -0.5 [1 , 2] = -0.5 [3 , 2] = 0.5 [2 , 3] = -0.5 [4 , 3] = 0.5 [3 , 4] = -0.5 [5 , 4] = 0.5 [4 , 5] = -0.5 [6 , 5] = 0.5 [5 , 6] = -0.5 [7 , 6] = 0.5 ⋮ [194, 195] = -0.5 [196, 195] = 0.5 [195, 196] = -0.5 [197, 196] = 0.5 [196, 197] = -0.5 [198, 197] = 0.5 [197, 198] = -0.5 [199, 198] = 0.5 [198, 199] = -0.5 [200, 199] = 0.5 [1 , 200] = 0.5 [199, 200] = -0.5

1
2
3
4

v = ones(N-1)
D2 = 2 * I - diagm( 1 => v, -1 => v)
D2[1,end] = D2[end,1] = -1
D2 = sparse( (1/Δx) * D2 )

200×200 SparseMatrixCSC{Float64,Int64} with 600 stored entries: [1 , 1] = 200.0 [2 , 1] = -100.0 [200, 1] = -100.0 [1 , 2] = -100.0 [2 , 2] = 200.0 [3 , 2] = -100.0 [2 , 3] = -100.0 [3 , 3] = 200.0 [4 , 3] = -100.0 [3 , 4] = -100.0 [4 , 4] = 200.0 [5 , 4] = -100.0 ⋮ [196, 197] = -100.0 [197, 197] = 200.0 [198, 197] = -100.0 [197, 198] = -100.0 [198, 198] = 200.0 [199, 198] = -100.0 [198, 199] = -100.0 [199, 199] = 200.0 [200, 199] = -100.0 [1 , 200] = -100.0 [199, 200] = -100.0 [200, 200] = 200.0

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

1
2
3
4
5
6

function f(u)
    b1 = Φ \ (D2 * u)
    b2 = - Ψ(u)*u - ϵ^2 * (D1 * b1)
    v = Φ \ b2
    return v
end

さあ、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

RK (generic function with 1 method

注: 上の $f(u)$ と Runge-Kutta スキームを見るとわかるが，時間方向に 1ステップ計算をすすめるたびに $Φ x = b$ という形の連立一次方程式を 8回解くことになる． しかもこの係数行列 $Φ$ は定数行列なので，全体計算に入る前に一回だけ LU分解(口頭で説明する)をして $L, U$ を保存しておけば「原理的には」だいぶ計算量を節約できそうだ． 余裕のある人はチャレンジしてみよう．
       ちなみに，今回の場合は行列 $Φ$ のサイズが小さく，かつ，その形が綺麗なので, (なにもしないと毎回行われてしまう) LU分解計算にかかる時間は比較的小さく，(Juliabox 上で, 以下同様)benchmark をとると 一回あたり 200～250$\mu$s で済んでしまう， そのため，びっくりするほどの大きな効果は出ないような気がするが，計算時間は環境に大きく依存するため，実際にやってみないとわからないものだ． たとえば全体で 50000 ステップの計算(2 ソリトンの奴)では，せいぜい $8$ 回/step $\times 50000$ step $\times 2.5 \times 10^{-4}$ 秒 $= 100$秒程度の短縮になるはずだが， 実際やってみると $250$秒ほど短縮され，「思ったより大きな効果が出る」改善となった． なお，計算結果が短縮前のものとほぼ同じになることは確認してある．

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

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

using ProgressMeter

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

@showprogress for i in 1:10000
    global u, sq_u
    # 関数定義の中同様，for 文の中は変数は外と切り離されているのでやり取りするなら global 宣言を．
    
    u = RK(u)
    if i % 1000 == 0  # 1000ステップおきにデータを格納するとして
        sq_u = hcat(sq_u, u)
    end
end

Progress: 100%|█████████████████████████████████████████| Time: 0:01:22 ↑ sparse matrix を使わないとこの計算に 20分ぐらいかかるよ！ ↑ ちなみに教官の事務作業用 PC だと 20秒で終わるぞ(もちろん sparse を使って)．

さて、試しに結果の数字をちょろっと見てみよう．

1

sq_u

200×11 Array{Float64,2}: 1.0 0.956114 0.859635 0.757929 … -0.580117 -0.645216 -0.625907 0.999507 0.964243 0.871681 0.770853 -0.573451 -0.579576 -0.534052 0.998027 0.971626 0.883347 0.783617 -0.508412 -0.461452 -0.392802 0.995562 0.978236 0.894615 0.796119 -0.377172 -0.279482 -0.19012 0.992115 0.984049 0.905469 0.808248 -0.169058 -0.022292 0.0839696 0.987688 0.98904 0.915888 0.819914 … 0.124332 0.314359 0.429629 0.982287 0.993182 0.925853 0.831088 0.499019 0.715756 0.825764 0.975917 0.996448 0.935342 0.841832 0.924747 1.13569 1.22019 0.968583 0.99881 0.944332 0.852306 1.33456 1.49389 1.53251 0.960294 1.00024 0.952801 0.862742 1.6329 1.69736 1.67991 0.951057 1.00071 0.960723 0.873388 … 1.73175 1.6837 1.61794 0.940881 1.00019 0.968071 0.884423 1.59854 1.45788 1.3668 0.929776 0.998652 0.974816 0.895864 1.27639 1.08895 0.999009 ⋮ ⋱ ⋮ 0.929776 0.808423 0.690582 0.594236 0.769585 1.28384 1.12515 0.940881 0.823813 0.706128 0.608616 1.04664 1.28126 0.855324 0.951057 0.838707 0.721445 0.622893 … 1.22533 1.11959 0.53291 0.960294 0.853086 0.736524 0.63705 1.2575 0.842215 0.214245 0.968583 0.866932 0.751353 0.651065 1.13277 0.512982 -0.0643768 0.975917 0.880225 0.765922 0.664921 0.883991 0.189695 -0.287608 0.982287 0.892947 0.78022 0.678605 0.570071 -0.0907169 -0.454673 0.987688 0.905077 0.794236 0.692118 … 0.25013 -0.312712 -0.571913 0.992115 0.916595 0.807957 0.705475 -0.034018 -0.475303 -0.647469 0.995562 0.92748 0.821371 0.718705 -0.261407 -0.584205 -0.688191 0.998027 0.937709 0.834465 0.731843 -0.426528 -0.646296 -0.698245 0.999507 0.947262 0.847224 0.744916 -0.531544 -0.666295 -0.678427

プロットもしておこうか．

1

plot( sq_u[:,end])

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

1

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

確かに、解がソリトンの和に分解されているようだ． 以前から書いているように、KdV 方程式の解は複数のソリトンの重ね合わせ的なもので書けることがわかっていて、ソリトンはその「高さ」によって移動速度が異なるため(背の高いほうが速い)、時間発展に従って速いものが前に素早く移動し、遅いものはだんだんとりのこされる．いま計算した結果は、だんだん分離するソリトンを見ていることになる．

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

I1(u) = sum(u) * Δx

I2(u) = sum( u.^2 ) * Δx

function I3(u)

    r1 = sum( (u.^3)/3 )
    r2 = sum( diff(u).^2 )/(Δx^2)
    r2 += (u[1]-u[N])^2/(Δx^2) # 面倒だけど、こいつもいれておかないとね…

    return (r1 - ϵ^2 * r2)*Δx
end

早速、それぞれの量の時間発展を計算してみよう．

1
2

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

11-element Array{Float64,1}: -3.552713678800501e-17 4.218847493575595e-17 -3.885780586188048e-17 1.509903313490213e-16 2.5202062658991054e-16 3.930189507173054e-16 4.085620730620576e-16 4.522771046566732e-16 3.241851231905457e-16 8.659739592076221e-17 1.3322676295501878e-16

1

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

11-element Array{Float64,1}: 1.0
1.0000132741706318 1.0000733293628314 1.000410860119632 1.0027302184279063 1.00725470213154
1.0106929683065309 1.0123549890163022 1.0127862928285352 1.0125304701638247 1.0119191376981238

1

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

11-element Array{Float64,1}: -0.004776495659718537 -0.004776252782382773 -0.004773009508861906 -0.004629894977744229 -0.002258141021547786 0.002656221917552628 0.005634993867614888 0.006376597912978426 0.00593554630678348 0.005044334397350881 0.004019924444297658

有限要素法を用いた場合でも、とりあえず、この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

function OneSoliton(height, center, x) # 1-ソリトン解
    d = sqrt(height/12) / ϵ
    return height * sech(d*(x-center))^2
end

u0 = [ OneSoliton(1, 0.5, k*Δx) + OneSoliton(0.5, 1.2, k*Δx) for k in 0:N-1 ]

200-element Array{Float64,1}: 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 0.0001435825782857728 0.00018666564694146636 ⋮ 6.6198161347727465e-6 5.498642874933181e-6 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

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

1

plot(u0)

早速計算してみよう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

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

@showprogress for i in 1:50000
    global u, sq_u

    u = RK(u)
    if i % 1000 == 0
        sq_u = hcat(sq_u, u)
    end
end

Progress: 100%|█████████████████████████████████████████| Time: 0:06:47 ↑ これも，sparse matrix を使わないと 90分ぐらいかかる． さらに上の注意(Runge-Kutta 関数を作ったすぐあと)で述べたような改善をすると，2分40秒程度で済む． ↑ なお，教官の事務作業用 PC (Julia 1.3.1) だと改善前で 1分30秒, 改善後で 1分10秒ぐらいだ．

さて、ではプロットしてみよう

1

plot(sq_u, legend = :no)

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

1
2
3

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

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

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

この計算では保存量は？

1
2

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

51-element Array{Float64,1}: 0.26019771077247456 0.2601977107724746 0.26019771077247467 0.2601977107724747 0.2601977107724747 0.2601977107724748 0.2601977107724748 0.2601977107724747 0.26019771077247467 0.2601977107724748 0.2601977107724748 0.2601977107724748 0.26019771077247483 ⋮
0.2601977107724747 0.2601977107724747 0.2601977107724747 0.2601977107724747 0.2601977107724746 0.2601977107724746 0.2601977107724746 0.2601977107724746 0.2601977107724746 0.2601977107724746 0.2601977107724746 0.2601977107724746

1

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

51-element Array{Float64,1}: 0.1375433423024314 0.13754270920262435 0.137542286175933
0.13754201154911835 0.13754190373461697 0.1375417757311086 0.13754155693795617 0.13754131868182976 0.13754122301921037 0.1375412977631334 0.13754094828384766 0.1375408457338157 0.13754045492095127 ⋮
0.13746616515261365 0.13748166374751164 0.13749501052183474 0.13750572837471386 0.13751394706208606 0.13752076225368998 0.13752617403380427 0.1375303067704669 0.13753340033001385 0.1375354332605638 0.13753714805536496 0.13753859098567292

1

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

51-element Array{Float64,1}: 0.023946737178329753 0.02394606703409743 0.02394578640424663 0.023945615361873746 0.023945591951155375 0.023945565002799027 0.023945441068107042 0.02394525607950774 0.023945217235507823 0.02394539216534583 0.02394513671543305 0.023945221930367282 0.02394504733053443 ⋮
0.023905577959153854 0.023914469479960274 0.023922064491820806 0.023927921725972048 0.023932063894224254 0.023935580182224223 0.02393833582343165 0.02394036208340304 0.023941917946647535 0.023942730744415563 0.02394350204041324 0.02394426610819958

グラフでみてみよう．

1
2

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

1

plot( sq_I2 .- sq_I2[1] )

1

plot( sq_I3 .- sq_I3[1] )

これもこれまで同様、二つのソリトンが衝突するあたりで数値誤差が発生している様子だ．

まあ、おおよそ上手くいっていると言えよう．

Report No.11 有限要素法で計算

1. 4つ目の保存量チェック
   wikipedia の「保存量」項 を参照して4つ目の保存量 $I_4$ を調べて、これ(の近似値)が数値計算でどう変化するか、今回(有限要素法)の計算方法に対して調べてみよ．グラフも描こう．
   
   
2. 熱拡散方程式 \[ u_t = u_{xx} \] をディリクレ境界条件 ($u(0,t) = 0, u(L,t) = 100$) などの境界条件のもとで、有限要素法スキームを自分で構成して数値計算してみよう(やや難しい)．
   
3. 上の(Runge-Kutta 関数の定義すぐ後の)計算時間短縮にチャレンジしてみよう．sparse 行列に対して LU 分解を行う SparseArrays.lu 関数については，Julia の中で ?lu としてマニュアルを参照しよう．web のマニュアルは記述があまりに不足なので要注意だ．