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

grid Photo by the blowup on Unsplash

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

差分法の類は空間の次元が上がっていくと「どう離散化を定義するか」という問題に直面することになる(空間次元が 1次元だと実感しにくいが)． そこが学術的には大変面白いところなのだが，ユーザとしては少し荷が重いかもしれない．

そうした意 味で，次元やメッシュの歪みに強い、汎用性の高い方法として 有限要素法(FEM: finite element method) を扱ってみよう．

ちなみにこれは数学上のちょっとした「発明」に基づく解法で，そうそう思いつけるものではない，画期的な解法だ． そのつもりで学ぼう．

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

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

KdV 方程式を対象として

ここでもやはり、Koreweg de Vries 方程式として, 次の式

\begin{equation} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \epsilon^2 \frac{\partial^3 u}{\partial x^3} = 0 \end{equation}

を考えることにしよう． ただし $\epsilon$ はパラメータだ．

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

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

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

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

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

早速計算してみよう！

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

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

1const ϵ,L = 0.022, 2
2
3N = 200
4Δx = L/N
5Δt = 0.0001  
6# やってみるとわかるが、今回のスキームはちょっと繊細で
7# Δt を細かくしないとうまく動かない．

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

1u0 = [ 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

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

1using Plots
2plot(u0, marker = :circ)

svg

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

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

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

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

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

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

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

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

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

1Matrix{Float64}(I, 3, 3)

とか

1diagm(0 => ones(3))

とすると $3 \times 3$ の(実数値の)単位行列が作れるぞ． どちらかというと後者のほうがすっきりしているかな．

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

 1function p(u)
 2  n = length(u)
 3  v = similar(u)
 4    
 5  for i in 2:n-1
 6    @inbounds v[i] = ( u[i+1] + u[i] + u[i-1] ) * ( u[i+1] - u[i-1] ) / 6.0
 7  end
 8    
 9  v[1] = ( u[2] + u[1] + u[n] ) * ( u[2] - u[n] ) / 6.0
10  v[n] = ( u[1] + u[n] + u[n-1] ) * ( u[1] - u[n-1] ) / 6.0
11    
12  return v
13end

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

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

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

そして D2 行列の定義だ．

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

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

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

1function f(u)
2    b1 = Φ \ (D2 * u)
3    b2 = p(u) +  (tD1 * b1)
4    v = - Φ \ b2
5    return v
6end

さあ、Runge-Kutta 法を定義しよう．

1function RK(u)
2    f1 = f(u)
3    f2 = f(u + Δt/2 * f1)
4    f3 = f(u + Δt/2 * f2)
5    f4 = f(u + Δt * f3)
6    return u + Δt * (f1 + 2*f2 + 2*f3 + f4)/6
7end

LU 分解の結果を取っておいて再利用することで高速化
上の $f(u)$ と Runge-Kutta スキームを見るとわかるが，時間方向に 1ステップ計算をすすめるたびに $\Phi \mathbf{x} = \mathbf{b}$ という形の連立一次方程式を 8回解くことになる． しかもこの係数行列 $\Phi$ は定数行列なので，全体計算に入る前に一回だけ LU分解 (掃き出し法の途中で出てくる情報を保存することに相当)をして $L, U$ を保存しておくと良い2． そうすると，連立一次方程式を解く計算の大半の部分を再利用できるので，計算時間が大いに減るのだ．

具体的には，

1const F = lu(Φ)

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

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

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

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

とすると良い4． こうすると $x$ に $\Phi \mathbf{x} = \mathbf{b}$ の解が入っている，というわけだ． 余裕のある人はこの高速化に(レポートで)チャレンジしてみよう．

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

1using BenchmarkTools
2
3@benchmark x = Φ \ u0

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  44.000 μs …   4.140 ms  ┊ GC (min … max): 0.00% … 44.11%
 Time  (median):     47.300 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   60.138 μs ± 113.998 μs  ┊ GC (mean ± σ):  1.17% ±  0.66%

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

 Memory estimate: 86.91 KiB, allocs estimate: 55.

一回の計算におおよそ 50～60μs ぐらいかかるということがわかる．

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

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

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  11.400 μs …  2.642 ms  ┊ GC (min … max): 0.00% … 98.89%
 Time  (median):     12.100 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   15.541 μs ± 36.313 μs  ┊ GC (mean ± σ):  6.28% ±  3.19%

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

 Memory estimate: 53.48 KiB, allocs estimate: 53.

おおよそ 15μs なので，この方法を採用すれば 2回目以降の計算は普通の計算方法よりも 3倍以上速いということになる．

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

1# const F = lu(Φ)
2x = similar(u0)
3@benchmark ldiv!( x, F, u0 )

BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.244 μs …   7.056 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.311 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.344 μs ± 215.981 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

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

 Memory estimate: 0 bytes, allocs estimate: 0.

おおよそ 2.3μs なのでさらに6～7倍速い．普通の計算に比べると，2回目以降の計算は 20倍以上速いということになり，爆速と言って良い．

というわけで，同じ係数行列の連立一次方程式を複数回解くときは LU分解の結果を再利用すべき，と頭の片隅に入れておこう．

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

 1using ProgressMeter
 2
 3u = copy(u0)
 4sq_u = copy(u0)
 5
 6@showprogress for i in 1:10000
 7    u = RK(u)
 8
 9    # 1000ステップおきにデータを格納
10    if i % 1000 == 0  
11        sq_u = hcat(sq_u, u)
12    end
13end

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

ちなみに，さきに示した「高速化」を施した右辺の式 f が次のものだ．

 1function f(u)
 2#    b1 = Φ \ (D2 * u)
 3    b1 = similar(u)
 4    ldiv!( b1, F, D2 * u)
 5    b2 = p(u) +  (tD1 * b1)
 6#    v = - Φ \ b2
 7    v = similar( b2 )
 8    ldiv!( v, F, -b2 )
 9    return v
10end

この f を定義してからもう一度 using ProgressMeter で始まる先の計算をし直してみると，ProgressBar が 0:00:00 という表示で終わるぞ.
真面目に @benchmark で測ると，10000step の計算が 290ms ぐらいで終わっていて，37倍も速い！ ということが確認できる．

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

1sq_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

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

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

svg

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

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

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

\begin{equation} \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. \end{equation}

という感じで、これまで通りそれぞれプログラムする．

 1I1(u) = sum(u) * Δx
 2
 3I2(u) = sum( u.^2 ) * Δx
 4
 5function I3(u)
 6
 7    r1 = sum( (u.^3)/3 )
 8    r2 = sum( diff(u).^2 )/(Δx^2)
 9    r2 += (u[1]-u[N])^2/(Δx^2) # 面倒だけど、こいつもいれておかないとね…
10
11    return (r1 - ϵ^2 * r2)*Δx
12end

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

1n = size(sq_u)[2]
2sq_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

1sq_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

1sq_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-ソリトン解は次の式で書ける．

\begin{equation} u(x) = H sech^2( \sqrt{\frac{H}{12\epsilon^2}} (x-x_0) ) \end{equation}

そこで、少し離してこれを二つ配置することで 2-ソリトンを擬似的に初期配置する．

1function OneSoliton(height, center, x) # 1-ソリトン解
2    d = sqrt(height/12) / ϵ
3    return height * sech(d*(x-center))^2
4end
5
6u0 = [ 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

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

1plot(u0, mark = :circ)

svg

早速計算してみよう.

1u = copy(u0)
2sq_u = copy(u0)
3
4@showprogress for i in 1:50000
5    u = RK(u)
6    if i % 1000 == 0
7        sq_u = hcat(sq_u, u)
8    end
9end

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

ちなみにこれは先のように高速化した場合の計算時間だ．高速化していない元のプログラムだと 1分ぐらいかかるんじゃないかな．

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

1plot(sq_u, legend = :no)

svg

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

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

svg

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

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

この計算では保存量は？

1n = size(sq_u)[2]
2sq_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

1sq_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

1sq_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

グラフでみてみよう．

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

svg

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

1plot( sq_I2 .- sq_I2[1] )

svg

1plot( sq_I3 .- sq_I3[1] )

svg

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

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

レポート

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

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

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

  注意
　　近年はセキュリティ上の懸念から，実行形式のプログラムなどをメールに添付するとそのメールそのものの受信を受信側サーバが拒絶したりする． そういうことを避けるため，レポートをメールで提出するときは添付ファイルにそういった懸念のあるファイルが無いようにしよう．

まあ要するに，レポートは pdf ファイルにして，それをメールに添付して送るのが良い ということだと思っておこう．

  1. 4つ目の保存量チェック
    wikipedia の「保存量」項 を参照して4つ目の保存量 $I_4$ を調べて、これ(の近似値)が数値計算でどう変化するか、今回(有限要素法)の計算方法に対して調べてみよ．グラフも描こう．

  2. 熱拡散方程式 \begin{equation} u_t = u_{xx} \end{equation} をディリクレ境界条件 ($u(0,t) = 0, u(L,t) = 100$) などの境界条件のもとで、有限要素法スキームを自分で構成して数値計算してみよう. 境界条件を有限要素法で実現するかがおそらく少し難しいだろう． 書籍などの資料等にあたるもヨシ，自分で考えるのもヨシだ．
    ヒントとしては，「境界に近いところだけ基底関数を特別なものにする」という感じかな．

  1. 素直に $A^{-1}$ と書けばよいではないかと思うかもしれないが，特殊な場合以外は数値計算では逆行列は用いない．というか，(サイズが小さい場合はともかく普通は)逆行列を用いたくても計算量的に無理なのだ．そういう意味で，気軽に逆行列を書いてしまう純粋数学の「癖」は罪深い… ↩︎

  2. LU 分解して $L$, $U$ を保存しておくことは，実質的に $A^{-1}$ を保存しておくことに相当する． ↩︎

  3. L, U, p, q, Rs というパラメータはそれぞれ以下のようなものである．まず L, U は LU 分解した行列で，p は Ax = b の b の順番を入れ替えた添字で，q は x の順番を入れ替えた添字．そして Rs は解の要素それぞれの縮小率に相当する．詳しくは ?lu としてマニュアルを見よう． ↩︎

  4. 理論上は速いんだけどきちんと実装しないと遅い，という話はたくさんある．こういう点に抜かり無いのが Julia の嬉しいところだ． ↩︎