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

有限要素法(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 法などを使えばよいだけだ．

早速計算してみよう！

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


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.995562 0.992115 0.987688 0.982287 0.975917 0.968583 0.960294 0.951057 0.940881 0.929776 ⋮
0.929776 0.940881 0.951057 0.960294 0.968583 0.975917 0.982287 0.987688 0.992115 0.995562 0.998027 0.999507

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

1
2
3


using Plots
gr()
plot(u0)

svg

スキームに出現する行列の類を作っていく

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)

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

 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

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

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

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[:,11])

svg

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

1

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

svg

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

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

これまで同様、KdV 方程式の保存量をチェックしよう．単純な方から再度挙げてみると、

$$ \left\{ \begin{array}{rcl} \displaystyle I_1 & := & \int_0^L u dx, \cr \displaystyle I_2 & := & \int_0^L u^2 dx, \cr \displaystyle I_3 & := & \int_0^L \left\{ u^3 / 3 - \epsilon^2 (u_x)^2 \right\}dx, \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

sq_I1 = [ I1(sq_u[:,k]) for k in 1:11 ]

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:11 ]

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:11 ]

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)

svg

早速計算してみよう.

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)

svg

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

1
2
3


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

svg

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

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

この計算では保存量は？

1

sq_I1 = [ I1(sq_u[:,k]) for k in 1:51 ]

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:51 ]

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:51 ]

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

plot(sq_I1 .- 0.260198 , marker = :circle)

svg

1

plot( sq_I2 .- 0.137543, marker = :circle)

svg

1

plot( sq_I3 .- 0.0239467, marker = :circle)

svg

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

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

Report No.19 自分で有限要素法(ややチャレンジング)

熱拡散方程式

$u_t = u_{xx}$

をディリクレ境界条件 ($u(0,t) = 0, u(L,t) = 100$) などの境界条件のもとで、有限要素法スキームを自分で構成して数値計算してみよう．