Photo by the blowup on Unsplash

有限要素法(FEM: Finite Element Method)とは

差分法の類は空間の次元が上がっていくと「どう離散化を定義するか」という問題に直面することになる(空間次元が 1次元だと実感しにくいが)． そこが学術的には大変面白いところなのだが，ユーザとしては少し荷が重いかもしれない． そうした意味で，次元やメッシュの歪みに強い、汎用性の高い方法として 有限要素法(FEM: 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 \]

を考えることにしよう．

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

詳細は上のノートに記載してあるが、ここでは解を空間方向には区分的一次多項式(つまり折れ線グラフだ)で近似することにして、それぞれのポイントの高さを集めたベクトルを $\boldsymbol{u}$ として、次のような近似式を得ている． 時間方向はそのままのため，これは一種の method of line の形である．

\[ \frac{d \boldsymbol{u}}{dt} = - \Phi \, \backslash \, \left\{ \, \boldsymbol{p}(\boldsymbol{u}) + \tilde{D}_1 \left( \, \Phi \,\backslash\, ( D_2 \boldsymbol{u} ) \, \right) \right\} \]

ただし、$\Phi, \boldsymbol{p}(\boldsymbol{u}), \tilde{D}_1, D_2$ はそれぞれ上のノート中で定義されている行列やベクトル値関数で，また， $A\, \backslash\, b$ は連立一次方程式 $A \boldsymbol{x} = b$ を解いて $\boldsymbol{x}$ を求めることを意味している¹．

よって，あとは Runge-Kutta 法などを使えば計算できてしまう． 今回はそうしよう．

早速計算してみよう！

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

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

1
2
3
4
5
6
7

const ϵ,L = 0.022, 2

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

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

1

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

200-element Array{Float64,1}:
 0.9995065603657316
 0.9980267284282716
 0.99556196460308
 0.9921147013144779
 0.9876883405951378
 ⋮
 0.9921147013144778
 0.99556196460308
 0.9980267284282716
 0.9995065603657316
 1.0

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

1
2

using Plots
plot(u0)

ふむ、いつもの初期値だな．ということで確認がとれたので次へ進もう．

次に，線形代数計算と，疎行列(「すかすか = ゼロ値な要素が多い」行列のこと)を扱うパッケージの利用宣言をしておこう． なお，この2つのパッケージは標準でインストールされているのでインストール作業は不要だ．

ここで出てくる SparseArrays package は疎な行列やベクトルを扱うのに適したパッケージで，適した問題にはメモリも計算量も少なくて済むことが多いので積極的に使っていこう．

1
2
3
4
5

using LinearAlgebra
# 線形代数の少し突っ込んだ機能はこのパッケージで．

using SparseArrays
# すかすかな行列やベクトルを扱うときはこれ．

さて，次にスキーム中に出現する行列を順番に作っていこう．最初は Φ 行列からだ．

1
2
3
4
5
6

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
⋮
[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$ の(実数値の)単位行列が作れるぞ．

さて話を戻して，次にベクトル関数 p を作るとしよう．

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

function p(u)
  n = length(u)
  v = similar(u)
    
  for i in 2:n-1
    @inbounds v[i] = ( u[i+1] + u[i] + u[i-1] ) * ( u[i+1] - u[i-1] ) / 6.0
  end
    
  v[1] = ( u[2] + u[1] + u[n] ) * ( u[2] - u[n] ) / 6.0
  v[n] = ( u[1] + u[n] + u[n-1] ) * ( u[1] - u[n-1] ) / 6.0
    
  return v
end

次に tD1 行列($\tilde{D}_1$のこと)だ．

1
2
3
4
5

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

200×200 SparseMatrixCSC{Float64,Int64} with 400 stored entries:
  [2  ,   1]  =  0.000242
  [200,   1]  =  -0.000242
  [1  ,   2]  =  -0.000242
  [3  ,   2]  =  0.000242
  [2  ,   3]  =  -0.000242
  ⋮
  [199, 198]  =  0.000242
  [198, 199]  =  -0.000242
  [200, 199]  =  0.000242
  [1  , 200]  =  0.000242
  [199, 200]  =  -0.000242

そして D2 行列の定義だ．

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
⋮
[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 = p(u) +  (tD1 * 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

LU 分解の結果を取っておいて再利用することで高速化
上の $f(u)$ と Runge-Kutta スキームを見るとわかるが，時間方向に 1ステップ計算をすすめるたびに $Φ x = b$ という形の連立一次方程式を 8回解くことになる． しかもこの係数行列 $Φ$ は定数行列なので，全体計算に入る前に一回だけ LU分解(口頭で説明する)をして $L, U$ を保存しておけば²「原理的には」だいぶ計算量を節約できそうだ．

具体的には，

1

const F = lu(Φ)

としてプログラム全体で一回だけ行列 Φ を LU 分解して結果を F に記憶させておく． そして，x = Φ \ b という連立一次方程式を解く箇所を 原理的には

1
2

x = similar(b)
x[F.q] = F.Rs .* ( F.U \ ( F.L \ b[F.p] ) )

に置き換えれば良い³．これで実質的に前進後退消去計算だけで連立一次方程式を解くことになる． ただ，これは内部的に計算に少し無駄がありちょっと遅いので，実際には LU 分解や QR 分解や cholesky 分解の結果を使って計算してくれる専用のコマンド ldiv! を使って

1
2

x = similar(b)
ldiv!( x, F, b )

とすると良い⁴． 余裕のある人はこの高速化に(レポートで)チャレンジしてみよう．

ちなみに，今回の場合 普通の PC で benchmark をとると Φ \ u0 は 40$\mu$s 程度で ldiv! 計算だとおおよそ 5$\mu$ なので，この工夫は「この部分のみ見るなら」 8倍程度の高速化効果をもつ． まあ，他の計算部分の時間もあるので，この高速化効果は「計算全体としては」もう少し緩いはず．あとは試してみるしかないね．

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

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

using ProgressMeter

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

@showprogress for i in 1:10000
    u = RK(u)

    # 1000ステップおきにデータを格納
    if i % 1000 == 0  
        sq_u = hcat(sq_u, u)
    end
end

 Progress: 100%|█████████████████████████| Time: 0:00:04 

さきに示した「高速化」を施してからすぐ上のこの計算をすると ProgressBar が 0:00:00 という表示で終わる(真面目に @benchmark で測ると，全部で 500ms ぐらいで計算が終わっている).

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

1

sq_u

200×11 Array{Float64,2}:
 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.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.0       0.956114  0.859635  0.757929     -0.580117  -0.645216   -0.625907

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

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}:
 -4.4408920985006264e-17
 -3.441691376337985e-17
  2.475797344914099e-16
  1.9317880628477724e-16
 -2.1094237467877975e-17
  1.9317880628477724e-16
  1.1546319456101628e-16
  2.9753977059954197e-16
  8.104628079763643e-17
 -8.43769498715119e-17
  2.3314683517128286e-17

1

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

11-element Array{Float64,1}:
 1.0
 1.000013274170632
 1.0000733293628319
 1.0004108601196324
 1.0027302184279068
 1.007254702131541
 1.0106929683065318
 1.0123549890163024
 1.012786292828537
 1.0125304701638254
 1.011919137698126

1

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

11-element Array{Float64,1}:
 -0.00477649565971854
 -0.004776252782382902
 -0.00477300950886177
 -0.00462989497774426
 -0.0022581410215480346
  0.0026562219175526635
  0.005634993867614746
  0.006376597912978284
  0.005935546306783834
  0.005044334397350881
  0.004019924444298155

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

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

@showprogress for i in 1:50000
    u = RK(u)
    if i % 1000 == 0
        sq_u = hcat(sq_u, u)
    end
end

Progress: 100%|███████████████████████████| Time: 0:00:20

これも高速化すると 0:00:02 ぐらいで終わるゾ．

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

1

plot(sq_u, legend = :no)

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

1
2
3

X = Δx:Δx:2       # 真面目に 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.26019771077247456
 0.26019771077247456
 0.2601977107724747
 ⋮
 0.26019771077247456
 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.13754228617593295
 0.1375420115491183
 0.13754190373461683
 ⋮
 0.13753030677046676
 0.13753340033001377
 0.1375354332605638
 0.13753714805536493
 0.1375385909856728

1

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

51-element Array{Float64,1}:
 0.023946737178329753
 0.02394606703409743
 0.02394578640424662
 0.023945615361873732
 0.023945591951155327
 ⋮
 0.023940362083402972
 0.023941917946647486
 0.023942730744415518
 0.023943502040413192
 0.0239442661081995

グラフでみてみよう．

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] )

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

まあ、十分に速いし，かなり上手くいっていると言えよう．

レポート

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

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

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

1. 4つ目の保存量チェック
   wikipedia の「保存量」項 を参照して4つ目の保存量 $I_4$ を調べて、これ(の近似値)が数値計算でどう変化するか、今回(有限要素法)の計算方法に対して調べてみよ．グラフも描こう．
   
   
2. 上の(Runge-Kutta 関数の定義すぐ後の)計算時間短縮にチャレンジしてみよう．どれだけ高速化できたか，@benchmark で計測しよう．
   
   
3. 熱拡散方程式 \[ u_t = u_{xx} \] をディリクレ境界条件 ($u(0,t) = 0, u(L,t) = 100$) などの境界条件のもとで、有限要素法スキームを自分で構成して数値計算してみよう. 境界条件を有限要素法で実現するかがおそらく少し難しいだろう． 書籍などの資料等にあたるもヨシ，自分で考えるのもヨシだ．
   ヒントとしては，「境界に近いところだけ基底関数を特別なものにする」という感じかな．