はじめに

偏微分方程式の数値解法として有限要素法は一般的な方法であり、多種多様な物理現象の計算に応用されています。

ただし、有限要素法が適用されるのは空間微分項に対してであり、時間微分項、すなわち時間発展（もしくは積分）の方法については有限差分法を用いることが（自分が知る限り）一般的です。複数の文献でも、時間積分については差分を導入することが特に注釈なく当たり前のように記述されています（例外もあるのですが、長くなるので後述します）。

有限要素法による時間積分の手法自体は30-40年前にはすでに提案され、数学者がまとめた有限要素法の書籍に短くはありますが記載されています¹,²。当時の計算技術では、この手法の高い計算コストが致命的な問題となり、多数派は差分法を採用した結果、今日の状況になったのかもしれません。

そこで、軽い気持ちで有限要素法での時間積分の方法について調べてみることにしました。有限要素法の初歩的な内容がわかっていれば理解できる程度までの記述とし、数学的に怪しい表現をしているので、気になる方は下に載せた文献を参照してください。

　弱形式の導出

まず、適当な例題を用意して、弱形式化された方程式を導出する。空間微分項の扱いは通常の有限要素法の作法で扱えるので、本記事では時間微分項の扱いのみを考慮する。その場合、時間に関する常微分方程式の解法を考えればよい。

問題設定

常微分方程式の典型的な例として、下記のスカラー $u$ の時間 $t$ に対する一次遅れ系を考える。

$\frac{\partial u}{\partial t}+Au=f, \text{ for } t\in[0,T], \text{ with } u(0)=v.$

ここで便宜的に適当な時間の最大値 $T$ を用意した。 $A$ については、ここではざっくりと正定値であることのみ定める。

　重み付き残差法

有限要素法の定石通り、まずは重み付き残差法を用いて上の式を弱形式で表現する。重み関数を $w, w(T)=0$ とし、

$\int_{0}^{T}{\left(\frac{\partial u}{\partial t}w+Auw-fw\right)}dt=0,$

が成り立つとする。このとき積分内の第一項は

$\int_{0}^T{\frac{\partial u}{\partial t}w dt}=\int_0^T{\left[\frac{\partial uw}{\partial t}-\frac{\partial w}{\partial t}u \right]dt}, \\=-[uw]_0^T-\int_0^T{\frac{\partial w}{\partial t}u dt}=-vw(0)-\int_0^T{\frac{\partial w}{\partial t}u dt}$

となり、基本境界条件が設定できる。以上をまとめると下記のようになる。

$\int_{0}^{T}{\left(-\frac{\partial w}{\partial t}u+Auw \right)}dt=vw(0)+\int_{0}^{T}{fw}dt.　\tag{I}$

離散化

有限要素法の通常の手順と同様、積分区間を有限の時間間隔で分割することを考え、分割区間内で変数 $u,w$ を多項式近似する。このとき $u,w$ に同じ近似関数を用いれば通常のGalerkin法となるが、全ての時間刻みの変数の値を自由度に持つことになり、本来は空間に対しても自由度を与えたいことを考えれば、自由度が膨大になりすぎて明らかに現実的ではない。解が時間に対して周期解になっていれば別だが、その場合は解を周波数空間に変換して解くほうが効率的であり、スペクトル要素法として整備されている。

不連続Galerkin法

上の問題を回避するための方法として、時間刻みの各点 $t_n$ において、前後の多項式が連続でないことを許容する（不連続Galerkin法）。

適当な時間 $t_n$ , $t_{n-1}$ を考え， $t_n-t_{n-1}=k_n$ として時間刻みを設定する。ここで $u,w$ ともに刻み幅の内部では連続とし、端点のみ不連続性を許容する。

すなわち、 $u^{n-1}$ と $u_+^{n-1}$ を用意し、下記の図のように定義する。

f:id:bigvalley:20200816215819p:plain

要するに

$\lim_{s \to 0+} u(t_{n-1}+s)=u_+^{n-1}$

を定義しておく。なお、Galerkin法の定義通り、 $w$ についても $u$ と同様の不連続性を与える。このとき、式(I)の左辺は

$\int_{0}^{T}{-\frac{\partial w}{\partial t}u}dt=\sum_{n=1}^{N-1}\left[-[uw]_{t_{n-1+0}}^{t_n}+\int_{t_{n-1}}^{t_n}{\frac{\partial u}{\partial t}w}dt \right],$

積分区間の間の不連続性を考慮し、

$[u]^n=u_{+}^{n}-u^n$

を導入すると

$\text{rhs}=\int_{0}^{T}{\frac{\partial u}{\partial t}w}dt+\sum_{n=1}^{N-1}\left([u]^nw \right)+(u_+^0 w^0),$

と表せる。式(I)に上の式を代入して整理すると、 $u$ $^n$ を除いて、各時間刻みで独立であることがわかる。そこで、適当な区間 $I_n=(t_{n-1},t_n$ ]について抜き出してみると、

$\int_{I_n}{\left(\frac{\partial u}{\partial t}w+Auw\right)}dt+\left([u]^{n-1} w_+^{n-1}\right)=\int_{I_n}{fw}dt,\tag{II}$

となる。最初の積分系の式の一部の積分区間のみを抜き出し、隣接区間との間の非独立成分を付け足した形になっている。隣接区間として右端を考慮していないが、未来（右端）から過去は変えられないと解釈すると、直感的に理解できるかもしれない。

多項式近似

最後に、 $u,w$ について区間 $I_n$ で多項式近似する。多項式の形はいろいろ考えられるが、後の計算の都合上、下記の $q-1$ 次の関数として表現する。

$u(t)=\sum_{j=0}^{q-1}\tilde{u}_j\left(\frac{t-t_{n-1}}{k_n}\right)^j.$

ここで $\tilde{u}_j$ は節点 $j$ に割り振られた値であり、括弧の項は通常の有限要素法ならば形状関数と呼ぶ（時間の多項式でも"形状"関数なのだろうか？）。0次近似ならば

$u(t)=\tilde{u}_0=u_+^{n-1}=u^n,$

1次近似ならば

$u(t)=\tilde{u}_0+\tilde{u}_1\left(\frac{t-t_{n-1}}{k_n}\right), u_+^{n-1}=\tilde{u}_0, \text{ }u^n=\tilde{u}_0+\tilde{u}_1.$

である．なお $w$ についても同じ形の近似関数を定義する．

通常の有限要素法の手順ではあまりやらないが、計算の簡略化のため、先に $w$ のみ展開してみる。 $w_+^{n-1}=\tilde{w}_0$ に注意して、1次の場合について書き下すと

$\int_{I_n}{\left(\frac{\partial u}{\partial t}+Au\right)\left(\tilde{w}_0+\tilde{w}_1\frac{t-t_{n-1}}{k_n}\right)}dt+(u_+^{n-1}-u^{n-1})\tilde{w}_0\\ =\int_{I_n}{f\left(\tilde{w}_0+\tilde{w}_1\frac{t-t_{n-1}}{k_n}\right)}dt,$

となる。これが任意の $\tilde{w}_0,\tilde{w}_1$ で成立する条件は下記の通り。

$\int_{I_n}{\left(\frac{\partial u}{\partial t}+Au\right)}dt+(u_+^{n-1}-u^{n-1})=\int_{I_n}{f}dt \tag{III.1}$

$\int_{I_n}{\left(\frac{\partial u}{\partial t}+Au\right)\left(\frac{t-t_{n-1}}{k_n}\right)}dt=\int_{I_n}{f\left(\frac{t-t_{n-1}}{k_n}\right)}dt, \tag{III.2}$

なお、0次の場合は式(III.1)のみが成立する。2次以上を考える場合も同様の手順でいけそう。

離散式

0次近似、1次近似の場合について，離散式を導出してみる。

まず0次近似について、式(III.1)に上の多項式を代入して整理する。積分区間での $u$ の時間微分が0であることに注意すると，下記のようになる。

$Au^nk_n+(u^n-u^{n-1})=\int_{I_n}{f}dt.　 \tag{IV}$

左辺は陰的Euler法の形、右辺は区間内平均の形になっている。

1次近似の場合について、同様に式(III.1),(III.2)に多項式を代入すると、

$\tilde{u}_1+\tilde{u}_0+k_nA\tilde{u}_1+\frac{1}{2}k_nA\tilde{u}_1=u^{n-1}+\int_{I_n}{f}dt \tag{V.1}$

$\frac{1}{2}k_nA\tilde{u}_0+\frac{1}{2}k_n\tilde{u}_1+\frac{1}{3}k_nA\tilde{u}_1=\frac{1}{k_n}\int_{I_n}{(t-t_{n-1})f}dt, \tag{V.2}$

となる。外力項 $f$ は各自の問題設定に従って連続関数として与えればいいので省略する。以上で、最終的に時間発展に必要な離散式（代数方程式）を実装レベルまで書き下せた。

解析解との比較

有限要素法による時間積分を考え、積分区間内での多項式を0次または1次とした場合について検討した。この表現がどの程度適切か、解析解との比較によって考えてみる。 $f=0$ とすれば、もとの常微分方程式の解析解は、漸化式の形で

$u(t_n)=r(Ak_n)u(t_{n-1})$

ここで $r(x)=\exp(-x)$ である。この形で式(IV)を整理すると、

$r_0(x)=\frac{1}{1+x},$

式(V)について

$r_1(x)=\frac{1-\frac{x}{3}}{1+\frac{2}{3}x+\frac{1}{6}x^2},$

となる（この式は $\exp(-x)$ のPadé近似らしい[1])（未確認）。

せっかくなので差分法の場合も考える。 Euler陽解法の場合、

$\frac{u^n-u^{n-1}}{k_n}+Au^{n-1}=0$

だから、 $r_\text{Ex}(x)=(1-x)$ となる。

次に、Crank-Nicolson法を用いると、

$\frac{u^n-u^{n-1}}{k_n}+\frac{1}{2}A(u^n+u^{n-1})=0$

だから、

$r_c(x)=\frac{1-\frac{x}{2}}{1+\frac{x}{2}}$

となる。

近似関数のプロット

解析解とそれぞれの解法で得られた $r(x)$ を描画してみると下記の通り。 f:id:bigvalley:20200817093444p:plain Euler陽解法が解析解から大きく外れていくのはご愛敬。明らかに1次多項式の不連続Galerkin法の結果が最も解析解に近くなった。0次多項式での不連続Galerkin法と等価なEuler陰解法だと、予想通り数値粘性が強く出てきて、解析解より値の減少が少ない。 Crank-Nicolson法は安定かつ実装がシンプル、それでいて二次精度と、使い勝手がよい印象だったが、1次の不連続Galerkin法には勝てなかった。このことは $r(x)$ から予想がついていて、各 $r(x)$ について $x\rightarrow\infty$ とすると、解析解、陰的Euler, 不連続Galerkin法は0に漸近するのに対して、Crank-Nicolson法だと-1に収束する。熱流体分野ではよく使われている印象だったが、気を抜いて刻み幅を大きくすると発散方向に解の誤差が大きくなるということか。

まとめ

有限要素法による時間積分について調べてみました。1次の不連続Galerkin法でさえ、計算の自由度は差分による陰解法の2倍となりますが、低次の差分解法と比較して精度向上が見込めそうです。とはいえ、理屈はかなり明快なので、それなりに使われていてもよさそうなものですが…。

この方法の自然な拡張の一つに有限要素法による流体構造連成計算手法の代表格ともいえるspace-time法(c.f., ³)があり、こちらに研究の動向が移り、わざわざ時間項のみの扱いに着目する研究者や技術者がいなくなったのかもしれません。

数十年前に確立された手法を興味本位で調べて時間を費やしたことに少し後悔していなくもないですが、誰かの興味や好奇心に少しでも刺されば幸いです。

Galerkin Finite Element Methods for Parabolic Problems (Springer Series in Computational Mathematics (25))

作者:Thomee, Vidar
発売日: 2006/09/15
メディア: ハードカバー

Numerical Solution of Partial Differential Equations by the Finite Element Method (Dover Books on Mathematics)

作者:Johnson, Claes
発売日: 2009/01/15
メディア: ペーパーバック

Computational Fluid-Structure Interaction: Methods and Applications (Wiley Series in Computational Mechanics)

作者:Bazilevs, Yuri,Takizawa, Kenji,Tezduyar, Tayfun E.
発売日: 2013/02/11
メディア: ハードカバー

V. Thomée, Galerkin Finite Element Methods for Paraboloc Problems, Springer, 1984. ↩
C. Johnson, Numerical solutions of partial differential equations by the finite element method, Dover, 2009↩
Y. Bazilevs, T. Takizawa, and T. E. Tezduyar. Computational fluid‐structure interaction: methods and applications. Wiley, 2013↩

明後日の方向にただ全力

微分方程式の時間発展を有限要素法で解く（不連続Galerkin法）