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

有限要素法(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$ が得られるので、あとは Runge-Kutta 法などを使えばよいだけだ(空間方向だけ離散化した結果を使うと、先に紹介した method of line の一種になる)．

早速計算してみよう！

では早速計算してみよう．常微分方程式の解法には 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.999507 ⋮ 0.998027 0.999507

念のためにプロットして確かめよう

1
2

using Plots
gr()

としてから、juliabox で julia のバージョンが 0.6.4 の人は

1

default( dpi = 450 )

ともしておこう．(図の出力が小さくなってしまわない)他の人はこの 1行は不要だ．

さて、出力だ．

1

plot(u0)

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

次に、スキーム中に出現する行列のたぐいを作っていこう．

1
2
3
4

v = ones(N-1)
Φ = (2/3)*eye(N) + (1/6)*( diagm(v, 1) + diagm(v, -1))
Φ[1,N] = Φ[N,1] = 1/6
Φ = Δx * Φ

200×200 Array{Float64,2}: 0.00666667 0.00166667 0.0 … 0.0 0.0 0.00166667 0.00166667 0.00666667 0.00166667 0.0 0.0 0.0 0.0 0.00166667 0.00666667 0.0 0.0 0.0 0.0 0.0 0.00166667 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋮ ⋱ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.00166667 0.0 0.0 0.0 0.0 0.0 0.00666667 0.00166667 0.0 0.0 0.0 0.0 0.00166667 0.00666667 0.00166667 0.00166667 0.0 0.0 0.0 0.00166667 0.00666667

次に、中心差分を作ってから…

 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(vc, 0) + diagm(vf, 1) + diagm(vf, -1)
    M[1,end] = M[end,1] = corner

    return M
end

Ψ (generic function with 1 method)

1
2
3
4
5
6
7
8
9

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

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

200×200 Array{Float64,2}: 200.0 -100.0 0.0 0.0 0.0 … 0.0 0.0 0.0 -100.0 -100.0 200.0 -100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -100.0 200.0 -100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -100.0 200.0 -100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -100.0 200.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -100.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋮ ⋱ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … -100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 200.0 -100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -100.0 200.0 -100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -100.0 200.0 -100.0 -100.0 0.0 0.0 0.0 0.0 0.0 0.0 -100.0 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

f (generic function with 1 method)

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

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

…というところなのだが、以前から書いている通り今の juliabox はプログラムの実行速度が大変に遅くなっている． たぶん、以下の数字のままだと 3時間(！) ぐらいかかってしまうだろう．

そこで、juliabox を利用している人は以下のプログラムを「仮に」一回実行させてどれくらい時間がかかりそうかを見積もってからその実行をすぐに止めて(画面上の Run の右にある ■ を(何回か)押せば良い)、せいぜい 10分ぐらいで終わるように for 文の数字を調整して実行しよう(得られる結果はちょっとだけになってしまうが)． なおそのとき、そのあとの if 文の数字の調整も忘れないようにしよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

using ProgressMeter

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

@showprogress for i in 1:10000
    u = RK(u)
    if i % 1000 == 0  # 1000ステップおきにデータを格納するとして
        sq_u = hcat(sq_u, u)
    end
end

Progress: 100%|█████████████████████████████████████████| Time: 0:06:16 ← 普通の PC だとこれぐらいで終わる．

まあ、汎用性が高い解法な分、ちょっと計算に時間が掛かる格好だな．まあ、プログラムも速度を追求して書いていないので、こんなものか．

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

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}: -2.44249e-17 -1.33227e-16 -6.66134e-18 -2.25375e-16 -6.66134e-17 -5.66214e-17 -9.65894e-17 -1.85824e-16 -5.09592e-16 -6.85008e-16 -7.53841e-16

1

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

11-element Array{Float64,1}: 1.0 1.00001 1.00007 1.00041 1.00273 1.00725 1.01069 1.01235 1.01279 1.01253 1.01192

1

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

11-element Array{Float64,1}: -0.0047765 -0.00477625 -0.00477301 -0.00462989 -0.00225814 0.00265622 0.00563499 0.0063766 0.00593555 0.00504433 0.00401992

有限要素法を用いた場合でも、とりあえず、この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.00653e-6 1.04091e-5 1.35327e-5 1.75936e-5 2.28731e-5 2.97369e-5 3.86603e-5 5.02614e-5 6.53436e-5 8.49514e-5 0.000110443 0.000143583 0.000186666 ⋮ 6.61982e-6 5.49864e-6 4.56736e-6 3.7938e-6 3.15126e-6 2.61754e-6 2.17421e-6 1.80597e-6 1.5001e-6 1.24603e-6 1.035e-6 8.59703e-7

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

1

plot(u0)

早速計算してみよう.

ただし、juliabox での計算時間が云々…というのは先と同じだ．自分で調整しよう．

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:28:29

結構時間がかかるなあ．授業時はパラメータを自分で変更する必要があるだろうな．

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

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.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 ⋮ 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198 0.260198

1

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

51-element Array{Float64,1}: 0.137543 0.137543 0.137542 0.137542 0.137542 0.137542 0.137542 0.137541 0.137541 0.137541 0.137541 0.137541 0.13754 ⋮ 0.137466 0.137482 0.137495 0.137506 0.137514 0.137521 0.137526 0.13753 0.137533 0.137535 0.137537 0.137539

1

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

51-element Array{Float64,1}: 0.0239467 0.0239461 0.0239458 0.0239456 0.0239456 0.0239456 0.0239454 0.0239453 0.0239452 0.0239454 0.0239451 0.0239452 0.023945 ⋮ 0.0239056 0.0239145 0.0239221 0.0239279 0.0239321 0.0239356 0.0239383 0.0239404 0.0239419 0.0239427 0.0239435 0.0239443

グラフでみてみよう．

1
2

default( 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$) などの境界条件のもとで、有限要素法スキームを自分で構成して数値計算してみよう(やや難しい)．