grid 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 Vector{Float64}:
 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)

svg

無駄に思ってもこまめに確認を！
このように，当たり前のグラフを得るためのプロットは無駄な作業に思うかもしれない．しかし，いろいろ作業をしていると思わぬミスや間違いがあるものだ．それらに気づかずに先へ進むとかえって大きなムダやトラブルになりかねない．多少面倒でもこのように頻繁にいろいろ確認するように心がけよう．

さて，いつもの初期値だな．ということで確認がとれたので次へ進もう．

次に，前回も暑かったように，線形代数計算と，疎行列(「すかすか = ゼロ値な要素が多い」行列のこと)を扱うパッケージの利用宣言をしておこう．

1
2
3
4


using LinearAlgebra
# 線形代数の少し突っ込んだ機能はこのパッケージで．
using SparseArrays
# すかすかな行列やベクトルを扱うときはこれ．

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

1
2
3
4
5
6


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

200×200 SparseMatrixCSC{Float64, Int64} with 600 stored entries:
⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈
⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦

細かいことだが，行列の右上と左下の要素にも非ゼロな数字が入っているはず，ということが上の表示からもわかることに注意しておこう．

前回も出てきたが，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 = sparse( diagm( -1 => v, 1 => -v ) )
tD1[1,end] = 1/2
tD1[end,1] = -1/2
tD1 = ϵ^2 * tD1

200×200 SparseMatrixCSC{Float64, Int64} with 400 stored entries:
⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈
⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦

そして D2 行列の定義だ．

1
2
3
4


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

200×200 SparseMatrixCSC{Float64, Int64} with 600 stored entries: ⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈ ⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀ ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀ ⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦

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

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

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 )

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

ちなみに，今回の場合の benchmark を見ておこう．まず，工夫なしでそのまま計算する Φ \ u0 の計算時間は以下のような感じだ．

1
2
3


using BenchmarkTools

@benchmark x = Φ \ u0

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  42.900 μs …  2.969 ms  ┊ GC (min … max): 0.00% … 26.33%
 Time  (median):     45.100 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   48.788 μs ± 74.596 μs  ┊ GC (mean ± σ):  2.43% ±  1.61%

    ▂▃▄▆ ▁█▅▃▄                                                 
  ▂▄████▇█████▆▅▅▄▄▃▃▄▃▃▄▃▃▃▃▃▃▃▃▃▃▃▂▂▃▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂ ▃
  42.9 μs         Histogram: frequency by time        56.8 μs <

 Memory estimate: 86.91 KiB, allocs estimate: 55.

おおよそ 45μs ぐらいかかるということがわかる．

では次に，LU 分解してとっておいた結果を機械的に再利用した場合が以下のような感じだ．

1
2
3


const F = lu(Φ)
x = similar(u0)
@benchmark x[F.q] = F.Rs .* ( F.U \ ( F.L \ u0[F.p] ) )

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  12.100 μs …  1.091 ms  ┊ GC (min … max): 0.00% … 98.18%
 Time  (median):     13.000 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   14.929 μs ± 31.801 μs  ┊ GC (mean ± σ):  7.00% ±  3.25%

   ▄▇▇██▇▄▅▄▃▁▁         ▁ ▂▂▃▁▂▂▂ ▁   ▂▂▂ ▁▂▂ ▁▁              ▂
  ██████████████▇▅▅▅▅▄▅▇████████████████████████▇▆▇█▇▇▆▆▆▁▅▅▄ █
  12.1 μs      Histogram: log(frequency) by time      22.4 μs <

 Memory estimate: 53.17 KiB, allocs estimate: 30.

おおよそ 13μs なので，普通に計算するより3倍以上速いということになるな．

そして最後に，LU分解してとっておいた結果を高速に適用する場合だ．

1
2
3


const F = lu(Φ)
x = similar(u0)
@benchmark ldiv!( x, F, u0 )

BenchmarkTools.Trial: 10000 samples with 7 evaluations.
 Range (min … max):  4.614 μs … 90.143 μs  ┊ GC (min … max): 0.00% … 93.45%
 Time  (median):     4.729 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.905 μs ±  2.437 μs  ┊ GC (mean ± σ):  1.79% ±  3.44%

  ▃▄█▆▇▇▃ ▁▁▁▂▂▂▁                                            ▂
  █████████████████▆▅▆▆▆▅▄▅▅▄▄▃▄▁▅▅▅▄▃▃▃▄▅▃▄▅▄▄▅▅▁▄▄▁▁▁▃▅▄▄▄ █
  4.61 μs      Histogram: log(frequency) by time     6.77 μs <

 Memory estimate: 9.41 KiB, allocs estimate: 2.

おおよそ 4.7μs なのでさらに3倍ほど速い．LU分解していないケースに比べるとおおよそ 10倍速いので，爆速と言って良いな．というわけで，同じ係数行列の連立一次方程式を複数回解くときは LU分解の結果を再利用すべき，と頭の片隅に入れておこう．

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

 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 Matrix{Float64}:
 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])

svg

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

1

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

svg

確かに、解がソリトンの和に分解されているようだ．以前から書いているように、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 Vector{Float64}:
 -4.4408920985006264e-17
  2.9309887850104134e-16
  4.007905118896815e-16
  5.118128143521972e-16
  3.930189507173054e-16
  2.708944180085382e-16
  2.020605904817785e-16
  5.995204332975846e-16
  7.649436639667329e-16
  9.026113190202524e-16
  7.227551890309769e-16

1

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

11-element Vector{Float64}:
 1.0
 1.0000132741706318
 1.0000733293628314
 1.000410860119632
 1.0027302184279066
 1.0072547021315401
 1.0106929683065315
 1.012354989016303
 1.0127862928285378
 1.0125304701638278
 1.011919137698127

1

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

11-element Vector{Float64}:
 -0.00477649565971854
 -0.004776252782382674
 -0.004773009508861716
 -0.004629894977744063
 -0.002258141021547875
  0.002656221917552415
  0.005634993867615031
  0.006376597912978852
  0.005935546306784687
  0.005044334397352515
  0.004019924444299221

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

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:00:22

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

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

1

plot(sq_u, legend = :no)

svg

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

1
2
3


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

svg

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

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

この計算では保存量は？

1
2


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

51-element Vector{Float64}:
 0.26019771077247456
 0.26019771077247456
 0.2601977107724746
 0.2601977107724746
 0.2601977107724746
 ⋮
 0.2601977107724746
 0.26019771077247456
 0.2601977107724745
 0.2601977107724745
 0.26019771077247456

1

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

51-element Vector{Float64}:
 0.1375433423024314
 0.1375427092026244
 0.13754228617593303
 0.13754201154911827
 0.1375419037346168
 0.13754177573110835
 ⋮
 0.13753030677046696
 0.13753340033001388
 0.1375354332605638
 0.13753714805536493
 0.1375385909856729

1

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

51-element Vector{Float64}:
 0.02394673717832975
 0.023946067034097455
 0.023945786404246636
 0.02394561536187373
 0.02394559195115532
 ⋮
 0.023940362083403034
 0.023941917946647517
 0.02394273074441553
 0.023943502040413192
 0.023944266108199518

グラフでみてみよう．

1
2


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

svg

かなりよく保存されているようだ．

1

plot( sq_I2 .- sq_I2[1] )

svg

1

plot( sq_I3 .- sq_I3[1] )

svg

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

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

レポート

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

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

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

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

素直に $A^{-1}$ と書けばよいではないかと思うかもしれないが，特殊な場合以外は数値計算では逆行列は用いない．というか，(サイズが小さい場合はともかく普通は)逆行列を用いたくても計算量的に無理なのだ．そういう意味で，気軽に逆行列を書いてしまう純粋数学の「癖」は罪深い… ↩︎
LU 分解して $L$, $U$ を保存しておくことは，実質的に $A^{-1}$ を保存しておくことに相当する． ↩︎
L, U, p, q, Rs というパラメータはそれぞれ以下のようなものである．まず L, U は LU 分解した行列で，p は Ax = b の b の順番を入れ替えた添字で，q は x の順番を入れ替えた添字．そして Rs は解の要素それぞれの縮小率に相当する．詳しくは ?lu としてマニュアルを見よう． ↩︎
理論上は速いんだけどきちんと実装しないと遅い，という話はたくさんある．こういう点に抜かり無いのが Julia の嬉しいところだ． ↩︎

数値計算法基礎 (2021)

11. 有限要素法 for 偏微分方程式