Help us understand the problem. What is going on with this article?

スペクトル要素法事始:Chebyshev多項式を用いて一次元Helmholtz方程式を解く

スペクトル要素法:有限要素法とスペクトル法を通じた理解

TwitterではRunge-Kuttaに続いて有限要素法が流行っていますが,2018年の数値計算Advent calendarにもいくつか有限要素法を扱った記事がありました.

有限要素法とは 弱形式(積分形式)で定式化した支配方程式を領域を分割した「要素」ごとに基底関数で展開し,線形の行列方程式を構成して解く手法です.領域を要素に分割するため複雑な形状に対応でき,現実世界での様々な問題を解く際にきわめて強力です.しかしながら,基底関数が低次のため精度に限界があるという欠点があります.1
それに対し,スペクトル法と呼ばれる方程式に現れる物理量を高次の展開関数で近似して計算する手法は,有限要素法とある種相補的な性質を持ちます.つまり,領域全体を一つの展開関数で記述するため領域の形状や境界条件に制限が生じる代わりに,展開関数の次数を上げることで精度を飛躍的に向上させることができます.ごく簡単にまとめると,

  • 有限要素法は領域を要素に分割するためさまざまな領域形状で問題を解くことができるが,展開関数が低次のため(スペクトル法に比べると)精度が良くない
  • スペクトル法は高次の展開関数を用いるため高精度だが,領域全体で一つの展開関数を用いるため領域の形状などに対し制限を抱える

といえます.2

そこで,これらの手法を組み合わせて 複雑な領域にも対応できる一般性をもち,しかも精度の良い手法 が欲しくなります.この手法こそが今回紹介する スペクトル要素法 です.スペクトル法あるいは有限要素法をご存じの方にとっては,

  • 有限要素法から見ると,各要素の基底関数としてより高次(かつ性質の良いもの)を用いる
  • スペクトル法から見ると,領域を要素に分割した上で弱形式の支配方程式にスペクトル法を適用し,要素境界には弱形式による連続性を用いる

と書くと理解しやすいと思います.


この記事では,スペクトル要素法が初めて提唱された1984年の論文
A spectral element method for fluid dynamics: Laminar flow in a channel expansion
の1.3節に基づいて,Dirichlet境界条件のもとでの一次元Helmholtz方程式

\begin{align}
\left\{ \begin{array}{l}
  \frac{\mathrm{d}^2 u(x)}{\mathrm{d} x^2} -\lambda^2 u(x) = f(x), ~~~~~-1 \le x \le 1 \\
  u(-1) = u(1) = 0 \\
\end{array} \right.
\end{align}

を解くプログラムをJuliaで実装し,スペクトル要素法の一端を紹介したいと思います.3

スペクトル要素法のアルゴリズム

領域の要素分割

一次元領域 $[x_{\text{I}}, x_{\text{O}}]$ を $M$ の要素に分割し, $i$ 番目の要素を $N_{x}^{i}+1$ の離散点で補間することを考えます.
各要素が $x^i \in [a^i, b^i]$ の領域を担当するものとして,

\begin{align}
  [x_{\text{I}}, x_{\text{O}}] &= [a^1, b^1] \cup [a^2, b^2] \cup \cdots \cup [a^i, b^i] \cup \cdots \cup [a^M, b^M] \nonumber \\
  &= \bigcup_{i=1}^{M} [a^i, b^i], ~~~~~~~~(a^1 = x_{\text{I}}, b^M = x_{\text{O}})
\end{align}



のようにします.


いま,ある要素を「大域座標」 $x \in [a^i, b^i]$ から「局所座標」 $\bar{x}^i \in [-1,1]$ に写すような写像を考えられます.
この局所座標を $\cos$ を用いて $N_{x}^i+1$ 点に離散化すると,

  \bar{x}_{j}^i = \cos \frac{\pi j}{N_{x}^i}, ~~~(j = 0, 1, \dots N_{x}^i)

のように書けます.


局所座標における離散点 $\bar{x}^i_j \in [-1,1]$ と,対応する大域座標 $x_q \in [a^i, b^i]$ には 要素境界のオーバーラップを除いて 一対一の関係があります.

Chebyshev多項式と物理量の近似

Chebyshev多項式は,

T_n(\theta) = \cos (n\theta)

のように表される関数で,基底関数としてよい性質を持っています.一次から五次までのChebyshev多項式を図示すると,


のようになります.Chebyshev多項式についてまとめられている資料をいくつか挙げておきます.

さて, $i$ 番要素での物理量 $u^i(\bar{x})$ をChebyshev多項式 $T_n$ を使って表します.$i$ 番要素の $j$ 次展開係数を $u_j^i$ と書くことにして,

\begin{align}
u^i(\bar{x}) &= \sum_{j=0}^{N_x^i} u_j^i h_j^i(\bar{x}^i) \\
h_j^i(\bar{x}^i) &= \frac{2}{N_i^x} \sum_{j=0}^{N_x^i} \frac{1}{\bar{c}_j \bar{c}_n} T_n(\bar{x}_j^i) T_n(\bar{x}^i)
\end{align}

\begin{align}
\bar{c}_k =\left\{ \begin{array}{ll}
  1 & (k\neq0, N_{x}^i) \\
  2 & (k=0, N_{x}^i) \\
\end{array} \right.
\end{align}

となります.$h_j^i$ は $i$ 番要素以外でゼロであり,Lagrange補間の性質

h_j^i(\bar{x}_k^i) = \delta_{jk}

を満たしています.このことから, 展開関数であるChebyshev多項式の次数=要素内の離散点数 であることが分かります.

次節以降で,Helmholtz方程式(の弱形式による定式化)をこの多項式によって展開し,解くべき一次方程式を定めます.

Helmholtz方程式の定式化

弱形式による定式化

Helmholtz方程式をスペクトル要素法で近似する(つまり,有限要素法と同じく弱形式で定式化し弱解を求める)ことは,次の言明と同値です.

\delta I = 0, ~~ I = \int_{x_{\text{I}}}^{x_{\text{O}}} \Bigl[ -\frac{1}{2} \Big( \frac{\mathrm{d} u}{\mathrm{d} x} \Big)^2 - \frac{\lambda^2}{2} u^2 -uf \Bigr] ~\mathrm{d}x

を満たす$u(x)$を見つける.

これを正当化する議論については,再掲ですが @t_kemmochi さんによる記事や書籍等を参照してください.

要素方程式

Helmholtz方程式の弱形式に対する要素方程式(The elemental equations)を次式で定めます.

\begin{align}
\int_{-1}^1 \Bigl[ -\frac{1}{2} \Big( \frac{\mathrm{d} u}{\mathrm{d} x} \Big)^2 - \frac{\lambda^2}{2} u^2 \Bigr] ~\mathrm{d}x &= \int_{-1}^1 uf  ~\mathrm{d}x \\
( A_{jk}^i -\lambda^2 B_{jk}^i ) u_{k}^i &= B_{jk}^i f_{k}^i
\end{align}

下付き添字のみ総和規約を用いることに注意します.

要素行列は以下より定まります.

\begin{align}
  A_{jk}^i &= \frac{-2}{L^i (N_{x}^i)^2} \frac{4}{\bar{c}_{j}\bar{c}_{k}}
  \sum_{n=0}^{N_{x}^i} \sum_{m=0}^{N_{x}^i} \frac{1}{\bar{c}_{n}\bar{c}_{m}} T_{n}(\bar{x}_j^i) T_{m}(\bar{x}_k^i) a_{nm}; \\
 a_{nm} &= \int_{-1}^{1} \Bigl\{ \frac{\mathrm{d} T_{n}}{\mathrm{d} x} \frac{\mathrm{d} T_{m}}{\mathrm{d} x} \Bigr\} ~\mathrm{d}x, \\
 B_{jk}^i &= \frac{L^i}{2 (N_{x}^i)^2} \frac{4}{\bar{c}_{j}\bar{c}_{k}} \sum_{n=0}^{N_{x}^i} \sum_{m=0}^{N_{x}^i} \frac{1}{\bar{c}_{n}\bar{c}_{m}} T_{n}(\bar{x}_j^i) T_{m}(\bar{x}_k^i) b_{nm}; \\
 b_{nm} &= \int_{-1}^{1} \{ T_n T_m \} ~\mathrm{d}x \nonumber
\end{align}

ここでChebyshev多項式の性質から,

\begin{align}
  a_{nm} &= \int_{-1}^{1} \Bigl\{ \frac{\mathrm{d} T_{n}}{\mathrm{d} x} \frac{\mathrm{d} T_{m}}{\mathrm{d} x} \Bigr\} ~\mathrm{d}x \\
  &= 0, &~~n+m ~~\mathrm{odd} \\
  &= \frac{nm}{2} [J_{|(n-m/2)}-J_{|(n+m/2)}], &~~n+m ~~\mathrm{even} \\
  J_{k} &= -4 \sum_{p=1}^k \frac{1}{2p-1}, &k \ge 1, \\
    &= 0, &k = 0, \\
 b_{nm} &= \int_{-1}^{1} \{ T_n T_m \} ~\mathrm{d}x \nonumber \\
  &= 0, &~~n+m ~~\mathrm{odd}\\
  &= \frac{1}{1-(n+m)^2} +\frac{1}{1-(n-m)^2}, &~~n+m ~~\mathrm{even}
\end{align}

です.これらの性質については先ほど挙げた竹広先生のノート吉田先生のノートを参照してください.
一見とても複雑な式の羅列ですが,前に述べたChebyshev多項式による物理量の近似式と見比べると,弱形式をChebyshev展開で定義したものであることが分かると思います.

大域方程式

各要素で定義された要素方程式を直接剛性法(Direct stiffness method)によって足し上げると,系全体の方程式,

  C_{pq} u_{q} = B_{pq} f_{q}

を得ることができます.直接剛性法における要素行列の足し上げを $B_{jk}^i$ から $B_{pq}$ への演算を例に見てみます.この操作の模式図を示すと,


のようになります.図から分かるように,同じ大域座標に対応する要素境界を意味する$B_{N_{x}^i k}^i$行と$B_{0 k}^{i+1}$行,$B_{j N_{x}^i}^i$列と$B_{k 0}^{i+1}$列が重ね合わされています.この処理は有限要素法においても全く一緒ですが,要素内部に相当する行列要素が(重ね合わされず)そのまま残っている点が異なります.
結局,$(N_{x}^i +1) \times (N_{x}^i +1)$の大きさの要素行列 $B_{jk}^i$ を $M$ 要素分足しあげたとき,大域行列 $B_{pq}$ の大きさは $\big(\sum_{i=1}^{M}N_{x}^i +1) \times (\sum_{i=1}^{M}N_{x}^i +1\big)$ となります.この行列の大きさが大域離散点数と一致していることが分かると思います.

境界条件

今回の実装において,Dirichlet境界条件は大域係数行列の対応する行を0とする行列圧縮(Matrix condensation)によって実現できます.この様子を模式図で示すと,


のようになります.

また,ゼロNeumann境界条件は係数行列をそのまま用いることで自然に導入(Naturally imposed)されることが弱形式の定式化から分かります.

実装と計算デモ

今回は,原論文でテスト問題として使われていたのと同じ条件で計算します.4

名前 条件
領域 $x \in [-1,1]$
定数項 $\lambda^2 = 0$
外力項 $f = \cos(\pi x+\pi /4)$
上流・下流における境界条件 $u[-1] = 0, u[1]=0$

このときの問題と解析解を示すと,

\begin{align}
\left\{ \begin{array}{l}
  \frac{\mathrm{d}^2 u(x)}{\mathrm{d} x^2} = \cos(\pi x+\pi /4), ~~~~~-1 \le x \le 1 \\
  u(-1) = u(1) = 0 \\
\end{array} \right.
\end{align}
u(x) = \frac{\sin(\pi x) - \cos(\pi x) - 1}{\sqrt{2} \pi^2}

となります.

実装はgithubに置いてあります.
Helmholtz_SEM_Julia
メインの処理は以下の通りです.

## Main
set_initial_condition(param_,var_)
  # 計算に利用する配列の初期化
set_global_x(param_,var_)
  # 大域座標と局所座標の間の変換を定義
set_B_local(param_,var_)
set_C_local(param_,var_)
  # 要素行列を定義
set_B_global(param_,var_)
set_C_global(param_,var_)
  # 大域行列を定義(直接剛性法による足しあげ)
set_f_local(param_,var_)
  # 局所座標での外力項を定義
set_f_global(param_,var_)
  # 大域座標での外力項を定義(直接剛性法による足しあげ)
set_RHS(param_,var_)
  # 大域右辺ベクトルを定義
var_.u = solve_LineraEquation(param_,var_)
  # 線形方程式を反復解法で解き,解を求める

結果

様々な要素数 $M$ と要素内離散点数 $N_x^i$ に対する数値解を,解析解と共に示します.

  • 要素数 $M=2$ ,要素内離散点数(=Chebyshev多項式の次数) $N_x^i=3$

  • 要素数 $M=4$ ,要素内離散点数 $N_x^i=6$

  • 要素数 $M=8$ ,要素内離散点数 $N_x^i=10$

計算点数( $\leftrightarrow$ 要素数,要素内離散点数)を増加させると,数値解が解析解とよく一致していくことがわかります.

要素数/離散点数による収束

解に対する要素数および要素内離散点数の影響をみていきます.

  • 要素数 $M=2,3,4$ ,要素内離散点数 $N_x^i=6$ に固定

  • 要素数 $M=2$ に固定,要素内離散点数 $N_x^i=3,5,7,9$

離散点の増加に伴い,それ以外の補間誤差(interpolant error)が減少していっていることが分かります.また,原論文では要素数を固定したとき要素内離散点数(=Chebyshev多項式の次数)に対して $L^\infty$ 誤差が指数的に減少することが述べられています.5

スペクトル要素法:応用

ここまで,スペクトル要素法が初めて提唱された論文に基づいて実装をおこなってきました.そのため,現在実際に使われている手法は今回実装したものとはかなり異なります.6
最近のスペクトル要素法について勉強したいという方は,例えば,
Spectral/hp Element Methods for Computational Fluid Dynamics
を参照してください.さらに最新の情報としては,この書籍の著者の一人でもあるS. Sherwin教授らのグループによるレビュー論文,
Spectral/hp element methods: recent developments, applications, and perspectives(Arxiv)
があります.

また,流体解析ソルバとしてスペクトル要素法を用いているものに,

  • Nektar++
    • University of UtahのM. Kirby教授,Imperial College LondonのS. Sherwin教授らの主導するプロジェクト
  • semtex
    • Monash UniversityのH. Blackburn教授によるプロジェクト.
    • 3次元Cartesian座標系に加え,Fourier展開を用いた一次元周期性のある系にも対応している
  • NEK5000

などがあります.ぜひ試してみて(そしてインストール方法や使い方をQiitaにまとめてみて)ください!

参考文献

Anthony Patera, A spectral element method for fluid dynamics: Laminar flow in a channel expansion
G. Karniadakis and S. Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics
竹広真一,チェビシェフ関数変換法
吉田茂生,直交多項式
チェビシェフ多項式による関数近似:Juliaで眺めてみる。
Juliaで数値計算 その1:コードサンプル〜基本的計算編〜
Julia言語と Plots + GR で複素関数のgifアニメーションを作る


  1. 要素数を増やした場合の誤差が代数的に減少する,という方が適切かもしれません. 

  2. これはごく大雑把な「お気持ち」程度の定義です.厳密な定義は書籍・論文等にあたってください.また,私が誤解している点がありましたらご指摘ください. 

  3. なお,執筆者の力不足から関数解析を用いたスペクトル要素法の基礎づけについては満足な議論ができておりません.申し訳ありません. 

  4. $\lambda^2=0$ としているためPoisson方程式になっており記事タイトルから外れていますが,実装はHelmholtz方程式としておこなっています 

  5. ここまでこの記事で示したかったのですが,原論文に基づいて $L^\infty$ を計算してみても論文通りの図が得られませんでした.バグ,あるいは私自身の理解が間違っている可能性があるため,この部分は宿題としたいと思います(いつかきちんと完成させたいです). 

  6. 例えば,展開関数としてはChebyshev多項式ではなくJacobi多項式やLegendre多項式(それらをさらに拡張したもの)などが使われますし,ある方向に周期性をもつ系に対してはFourier展開と2次元スペクトル要素法を組み合わせて解く,といったこともおこなわれています. 

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away