1次元領域に対する有限要素法
この記事では、閉区間\([0, 1]\)上の境界条件 \(u(0) = u(1) = 0\) 付きのラプラシアンの固有値問題 \(-u''(x) = \lambda u(x)\) を有限要素法を用いて数値的に解く。この問題の解は三角関数を用いて簡単に表せることがよく知られているので有限要素法で数値計算する必要は全く無いが、有限要素法を理解するために手を動かすのにちょうどいいやさしい問題なので、ここに計算ノートを記しておく。
なお、MFEMをつかって2次元領域についての問題を解く手順については、前回の記事で書いた。
ラプラシアンの固有値問題
まずは、一般の\(d\)次元でのラプラシアンの固有値問題を考える。すなわち、\(d\)次元有界領域\(\Omega \subset \mathbb{R}^d\)に対して、\(\overline{\Omega}\)上の\(C^2\)級関数 \(u\colon \overline{\Omega} \to \mathbb{R}\)であって、
\[\begin{eqnarray} -\Delta u(x) &=& \lambda u(x) & (x \in \Omega), \\ u(x) &=& 0 & (x \in \partial \Omega) \end{eqnarray}\]を満たすものを考える。ここで、\(\Delta\)はラプラシアン\(\frac{\partial^2}{\partial x_1^2} + \cdots + \frac{\partial^2}{\partial x_d^2}\)。
この問題はたとえば、波動方程式
\[\begin{eqnarray} \frac{\partial^2 f}{\partial t^2}(x, t) &=& c^2\Delta f(x, t) & (x \in \Omega, t \in \mathbb{R}) \\ f(x, t) &=& 0 & (x \in \partial \Omega, t \in \mathbb{R}) \end{eqnarray}\]を変数分離すると現れる。
次元\(d=1\)の場合
領域の次元\(d\)が\(1\)の場合、すなわち\(u\colon [0,1] \to \mathbb{R}\)であって
\[\left \{ \begin{eqnarray} -u''(x) = \lambda u(x) &\quad& (0 < x < 1) \\ u(0) = u(1) = 0 & & \end{eqnarray} \right.\]を満たすものを求める問題の解はよく知られており、\(n \in \mathbb{N}\)に対して
\[\begin{eqnarray} \lambda &=& n^2\pi^2, \\ u(x) &=& \sin \left(n\pi x\right) \end{eqnarray}\]となる。
以下では、この\(d=1\)に制限した場合のみを扱う。
弱形式
有限要素法を適用するために、微分方程式を同値な別の形に変形する。
命題（ref. 『偏微分方程式の計算数理』命題3.2）
\(C^2\)級関数\(u\colon [0,1] \to \mathbb{R}\)で\(u(0)=u(1)=0\)であるものが\(-u''(x)=\lambda u(x)\)を満たすための必要十分条件は、\(v''(x) = \lambda v(x), v(0)=v(1)=0\)を満たす任意の\(C^2\)級関数\(v \colon [0,1]\to \mathbb{R}\)に対して、\(u(x)\)が
\[\int_0^1 u'(x)v'(x)dx = \lambda \int_0^1 u(x)v(x)dx\]を満たすことである。
証明
必要性：\(-u''(x)=\lambda u(x)\)の両辺に\(v(x)\)をかけて積分して\(-\int_0^1u''(x)v(x)dx = \lambda \int_0^1 u(x)v(x)dx\)を得る。部分積分と\(v(0)=v(1)=0\)の仮定から、左辺は\(-\int_0^1u''(x)v(x)dx = -\left[u'(x)v(x)\right]_0^1 + \int_0^1u'(x)v'(x)dx = \int_0^1u'(x)v'(x)dx\)と計算でき、必要性がわかる。
十分性：計算を逆にたどって、\(\int_0^1\left(u''(x)+\lambda u(x)\right)v(x)dx = 0\)が導ける。\(-u''(x) \neq \lambda u(x)\)と仮定して矛盾を導く。仮定より、ある点\(x_0\)で\(u''(x_0) + \lambda u(x_0) \neq 0\)となる。\(u\)は\(C^2\)級なので\(u''(x) + \lambda u(x)\)は連続であり、したがって\(x_0\)を含むある近傍上で\(u''(x) + \lambda u(x) > 0\)（または\(<0\)）となる。\(v\)としてその近傍内でのみ正の値をとる\(C^\infty\)級関数をとれば、\(\int_0^1\left(u''(x)+\lambda u(x)\right)v(x)dx > 0\)（または\(<0\)）となり、仮定に反する。（証明終わり）
以下では、微分方程式\(-u''(x)=\lambda u(x)\)の代わりに、この積分の方程式
\[\int_0^1 u'(x)v'(x)dx = \lambda \int_0^1 u(x)v(x)dx\]を使う。
有限要素法
計算機で計算可能にするために、閉区間\([0,1]\)を有限個の点で近似し、\(C^2\)級関数の集合 \(\{ v \colon [0,1] \to \mathbb{R} \mid v(0) = v(1) = 0\}\)を有限次元の線型空間で近似する。
まず、閉区間\([0,1]\)を線分の和に分割する。すなわち、\(n+1\)個の点\(P_0, P_1, \dots, P_{n-1}, P_n\)を\(0 = P_0 < P_1 < \cdots < P_{n-1} < P_n = 1\)を満たすようにとる。ここでは計算を簡単にするため、均等な分割 \(P_i = i / n \quad (0 \leq i \leq n)\)を使う。
境界である端点\(P_0, P_n\)を除く\(P_i \ (1 \leq i \leq n - 1)\)にそれぞれ対して、その点\(P_i\)でのみ\(1\)をとりその他の点\(P_j \ (j \neq i)\)で\(0\)をとる関数\(\varphi_i: [0, 1] \to \mathbb{R}\)を次のように定義する。
\[\varphi_i(x) = \left\{ \begin{eqnarray} nx - \left(i - 1\right) && \left(\frac {i - 1}{n} \leq x \leq \frac{i}{n}\right) \\ -nx + \left(i + 1\right) && \left(\frac{i}{n}\leq x\leq\frac{i+1}{n}\right) \\ 0 && \left(\text{otherwise}\right) \end{eqnarray} \right.\]
先ほど求めた積分の方程式を、有限集合 \(\varphi_1, \dots, \varphi_{n-1}\) が張る線型空間 \(\left\{a_1\varphi_1(x) + \cdots + a_{n-1}\varphi_{n-1}(x) \mid a_1, \dots, a_{n-1} \in \mathbb{R}\right\}\) 上に制限して考える。解\(u(x)\)を\(n-1\)個の実数\(u_1, \dots, u_{n-1}\)を用いて\(u(x) = \sum_{j=1}^{n-1} u_j\varphi_j(x)\)と表し、各\(i \ (1 \leq i \leq n - 1)\)に対して\(v(x) = \varphi_i(x)\)の場合に計算すると、
\[\begin{eqnarray} \int_0^1 \left(\sum_{j=1}^{n-1} u_j\varphi_j(x)\right)' \varphi_i'(x)dx = \lambda\int_0^1 \left(\sum_{j=1}^{n-1}u_j\varphi_j(x)\right)\varphi_i(x)dx \\ \end{eqnarray}\]となる。\(n-1\)次正方行列\(K\), \(M\)を、第\((i, j)\)成分がそれぞれ\(\int_0^1 \varphi_i'(x) \varphi_j'(x)dx\)、\(\int_0^1\varphi_i(x)\varphi_j(x)dx\)であるものとして定義する。\(n-1\)次元ベクトルを\(u = (u_1 \cdots u_{n-1})^\intercal\)と定義すれば、上の式は
\[Ku = \lambda Mu\]となる。いま、\(M\)は正定値対称行列（e.g. 『ラプラシアンの幾何と有限要素法』補題6.8）なので正則だから、上式は\(M^{-1}Ku = \lambda u\)と同値である。これは行列\(M^{-1}K\)の固有値および固有ベクトルを求める問題となっている。
ここで、行列\(K\)と\(M\)の値を具体的に計算しよう。
まず、\(\int_0^1\varphi_i'(x)\varphi_j'(x)dx\)は、\(\vert i - j\vert > 1\)のときには\(0\)となる。\(i=j\)のときは、
\[\begin{eqnarray} & &\int_0^1\varphi_i'(x)\varphi_i'(x)dx \\ &=& 2\int_0^{1/n} n^2dx \\ &=& 2n. \end{eqnarray}\]\(j = i + 1\)のときは、
\[\begin{eqnarray} & &\int_0^1\varphi_i'(x)\varphi_{i+1}'(x)dx \\ &=&\int_0^{1/n}(-n^2)dx \\ &=&-n. \end{eqnarray}\]
次に、\(\int_0^1\varphi_i(x)\varphi_j(x)dx\)についても、\(\vert i - j \vert > 1\)のときには\(0\)となる。\(i = j\)のときは、
\[\begin{eqnarray} &&\int_0^1\varphi_i(x)\varphi_i(x)dx \\ &=&2\int_0^{1/n}(nx)^2dx \\ &=&\frac{2}{3n}. \end{eqnarray}\]\(j = i + 1\)のときは、
\[\begin{eqnarray} &&\int_0^1\varphi_i(x)\varphi_{i+1}(x)dx \\ &=&\int_0^{1/n}nx(1-nx)dx \\ &=&\frac{1}{6n}. \end{eqnarray}\]数値計算
では、実際に有限要素法をPythonを使って実行してみよう。
まず、上で計算した通り、行列\(K, M\)をそれぞれ計算する。
import numpy as np
N_MESH = 1000
def calc_stiffness_matrix(n_mesh):
assert n_mesh >= 2
ret = np.zeros((n_mesh - 1, n_mesh - 1))
for i in range(n_mesh - 1):
ret[i][i] = 2 * n_mesh
for i in range(n_mesh - 2):
ret[i][i + 1] = -n_mesh
ret[i + 1][i] = -n_mesh
return ret
def calc_mass_matrix(n_mesh):
assert n_mesh >= 2
ret = np.zeros((n_mesh - 1, n_mesh - 1))
for i in range(n_mesh - 1):
ret[i][i] = 2 / (3 * n_mesh)
for i in range(n_mesh - 2):
ret[i][i + 1] = 1 / (6 * n_mesh)
ret[i + 1][i] = 1 / (6 * n_mesh)
return ret
stiffness_matrix = calc_stiffness_matrix(N_MESH)
mass_matrix = calc_mass_matrix(N_MESH)
次に、固有値と固有ベクトルを計算する。
先ほど行列\(M^{-1}K\)の固有値問題\(M^{-1}Ku = \lambda u\)を解けばいいと述べたが、数値計算の実用上は、逆行列\(M^{-1}\)の計算は避けて一般化固有値問題\(Ku = \lambda M u\)としたまま解くほうがいいと思われる。固有値・固有ベクトルの計算する関数は scipy.linalg.eigh や scipy.sparse.linalg.eigsh などがあり、残念ながら私は違いがよくわかっていないので適当に選んだ。
import scipy
values, vectors = scipy.linalg.eigh(stiffness_matrix, mass_matrix)
最初の3つの固有値を確認すると、以下のようになる。
print(values[:3])
[ 9.86961252 39.47854748 88.82709712]
計算された最初の3つの固有値が、実際の固有値\(1^2\pi^2\)、\(2^2\pi^2\)、\(3^2\pi^2\)に近い値をとっていることが見てとれる。
次に、最初の3つの固有関数をプロットすると以下のようになる。
import matplotlib.pyplot as plt
plt.plot(vectors[:, :3])
plt.savefig("solutions.png")
最初の3つの解\(\sin n\pi x\)、\(\sin 2n\pi x\)、\(\sin 3n\pi x\)に（定数倍をのぞいて）近い関数になっていることが観察できる。