第3章: 有限要素法の基礎 - 偏微分方程式の数値解法

学習目標

第1章の有限差分法(Finite Difference Method)、第2章のCrank-Nicolson法では、微分を格子点上の差分商で近似する立場を学びました。本章では、まったく異なる発想に立つ有限要素法(Finite Element Method、FEM)を学びます。有限要素法は、方程式を「弱形式(weak form)」に書き換え、領域を小さな要素に分割して、区分的な関数で解を近似します。この考え方は複雑な形状や不均一なメッシュに柔軟に対応でき、構造力学・伝熱・電磁場解析など幅広い分野で標準的な手法となっています。

本章を終えると、次のことができるようになります。

有限要素法と有限差分法の発想の違いを説明できる
1次元Poisson方程式の弱形式を、部分積分を用いて自分で導出できる
線形形状関数から要素剛性行列を構成し、全体行列へ組み立て(アセンブリ)できる
Dirichlet境界条件とNeumann境界条件の扱いの違いを理解する
NumPyだけを用いて1次元有限要素法ソルバーをゼロから実装し、厳密解と比較して収束次数を確認できる

読了時間の目安: 30〜35分 / 難易度: 中級(Intermediate) / コード例: 4 / 演習問題: 3

3.1 有限要素法とは

有限差分法は、微分方程式に現れる導関数を、格子点上の値の差分で直接置き換える方法でした。実装が直感的で、規則的な格子上では非常に強力です。一方で、複雑な境界形状や、局所的に細かい解像度が必要な問題では、規則格子の制約が弱点になります。

有限要素法(Finite Element Method、FEM)は、対象領域を 要素(element)と呼ぶ小区間(1次元では線分、2次元では三角形や四辺形)に分割し、各要素の端点である節点(node)での値を未知数として、解を区分的な基底関数(basis function)の重ね合わせで表現します。微分方程式そのものではなく、それを積分の形に緩めた弱形式(weak form)を満たすように未知数を決めるのが特徴です。

有限要素法の基本用語

要素(element): 領域を分割した小さな部分区間。1次元では隣り合う節点で挟まれた線分。
節点(node): 要素の端点。ここでの解の値が直接の未知数になります。
形状関数(shape function): 各要素内で解を補間する基底関数。本章では区分1次(線形)関数を用います。
剛性行列(stiffness matrix): 弱形式の左辺から生じる係数行列。構造力学に由来する呼び名です。

有限差分法との違いを整理すると次のようになります。

観点	有限差分法(FDM)	有限要素法(FEM)
近似の対象	各点での導関数	弱形式(積分方程式)
解の表現	格子点の値	形状関数の線形結合
メッシュ	規則格子が基本	不規則メッシュに柔軟
複雑形状	取り扱いが難しい	得意
境界条件	直接代入する	Dirichletは強制、Neumannは自然に入る

本章の題材: 説明を明確にするため、本章では1次元Poisson方程式 \( -u''(x) = f(x) \)、\( x \in (0,1) \)、境界条件 \( u(0)=u(1)=0 \) に絞って議論します。次元は低くても、弱形式・形状関数・アセンブリ・境界条件処理という有限要素法の骨格はすべてここに現れます。

3.2 弱形式と変分原理

出発点となる微分方程式そのものを強形式(strong form)と呼びます。

強形式

\[ -u''(x) = f(x), \quad x \in (0,1), \qquad u(0) = u(1) = 0. \]

弱形式を得るには、試験関数(test function) \( v(x) \) を導入します。 \( v \) は境界で \( v(0)=v(1)=0 \) を満たす十分滑らかな関数とします。強形式の両辺に \( v \) を掛け、区間 \( (0,1) \) で積分します。

\[ -\int_0^1 u''(x)\,v(x)\,dx = \int_0^1 f(x)\,v(x)\,dx. \]

左辺に部分積分(integration by parts)を適用します。

\[ -\int_0^1 u''\,v\,dx = -\bigl[\,u'\,v\,\bigr]_0^1 + \int_0^1 u'(x)\,v'(x)\,dx. \]

ここで境界項 \( [u'v]_0^1 \) は、試験関数が \( v(0)=v(1)=0 \) を満たすため消えます。したがって次の弱形式が得られます。

弱形式(変分形式)

すべての試験関数 \( v \)(\( v(0)=v(1)=0 \))に対して、

\[ \int_0^1 u'(x)\,v'(x)\,dx = \int_0^1 f(x)\,v(x)\,dx. \]

強形式では \( u \) に2階微分可能性を要求していましたが、弱形式では \( u \) と \( v \) の1階微分しか現れません。要求される滑らかさが「弱まる」ことが、弱形式という名前の由来です。この緩和のおかげで、区分1次関数のような折れ線でも解の候補として使えるようになります。

Galerkin法(Galerkin method)

未知関数を有限個の基底関数 \( \phi_j \) で \( u(x) \approx u_h(x) = \sum_{j} u_j\,\phi_j(x) \) と近似し、試験関数にも同じ基底 \( \phi_i \) を用いる方法をGalerkin法と呼びます。弱形式に代入すると、節点値 \( u_j \) に関する連立1次方程式 \( \sum_j \left(\int_0^1 \phi_i' \phi_j'\,dx\right) u_j = \int_0^1 f\,\phi_i\,dx \) が得られます。左辺の係数行列が剛性行列 \( K \)、右辺が荷重ベクトル \( F \) です。

3.3 形状関数と要素剛性行列

区間 \( [0,1] \) を \( n_e \) 個の要素に分割し、節点を \( 0 = x_0 < x_1 < \dots < x_{n_e} = 1 \) とします。各節点 \( x_j \) に対応する基底関数 \( \phi_j \) として、「\( x_j \) で1、隣の節点で0」となる三角形状の ハット関数(hat function)を用います。これが区分1次の形状関数です。

実装では、各要素を基準要素 \( \xi \in [-1, 1] \) に写して考えると便利です。基準要素上の線形形状関数は次の2つです。

\[ N_1(\xi) = \frac{1-\xi}{2}, \qquad N_2(\xi) = \frac{1+\xi}{2}. \]

\( N_1 \) は左端 \( \xi=-1 \) で1・右端で0、\( N_2 \) はその逆で、常に \( N_1+N_2=1 \)(1の分割)が成り立ちます。幅 \( h \) の要素上で、形状関数の導関数は \( \dfrac{d N_1}{dx} = -\dfrac{1}{h} \)、\( \dfrac{d N_2}{dx} = +\dfrac{1}{h} \) です。

弱形式の左辺を1要素分だけ取り出したものが要素剛性行列(element stiffness matrix)です。

要素剛性行列

\[ K^{e}_{ab} = \int_{x_L}^{x_R} \frac{dN_a}{dx}\,\frac{dN_b}{dx}\,dx = \frac{1}{h} \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix}. \]

導関数が要素内で一定なので、積分は幅 \( h \) を掛けるだけで求まります。

形状関数と要素剛性行列をコードで確認します。以下はすべて実行し、実際の出力を掲載しています。

コード例1: 形状関数と要素剛性行列

import numpy as np

def shape_functions(xi):
    """基準要素 [-1, 1] 上の線形形状関数 N1, N2 を返す"""
    N1 = (1.0 - xi) / 2.0
    N2 = (1.0 + xi) / 2.0
    return N1, N2

def element_stiffness(h):
    """幅 h の要素の要素剛性行列を返す"""
    return (1.0 / h) * np.array([[1.0, -1.0],
                                 [-1.0, 1.0]])

print("=== 形状関数の値 (1の分割 N1+N2=1 を確認) ===")
for xi in [-1.0, -0.5, 0.0, 0.5, 1.0]:
    N1, N2 = shape_functions(xi)
    print(f"xi={xi:+.1f}: N1={N1:.3f}, N2={N2:.3f}, sum={N1+N2:.3f}")

print()
print("=== 要素剛性行列 (h=0.25) ===")
print(element_stiffness(0.25))

=== 形状関数の値 (1の分割 N1+N2=1 を確認) === xi=-1.0: N1=1.000, N2=0.000, sum=1.000 xi=-0.5: N1=0.750, N2=0.250, sum=1.000 xi=+0.0: N1=0.500, N2=0.500, sum=1.000 xi=+0.5: N1=0.250, N2=0.750, sum=1.000 xi=+1.0: N1=0.000, N2=1.000, sum=1.000 === 要素剛性行列 (h=0.25) === [[ 4. -4.] [-4. 4.]]

どの評価点でも \( N_1+N_2=1 \) が成り立ち、形状関数が正しく1を分割していることが確認できます。 \( h=0.25 \) のとき \( 1/h = 4 \) なので、要素剛性行列の成分が \( \pm 4 \) になっている点も理論と一致します。

3.4 境界条件の処理と求解

各要素の要素剛性行列を、共有する節点の位置で足し合わせて全体の 剛性行列(global stiffness matrix) \( K \) を作る操作を アセンブリ(assembly)と呼びます。要素 \( e \)(節点 \( e \) と \( e+1 \) を持つ)の寄与を、全体行列の対応する行・列に加算していきます。隣り合う要素は節点を共有するため、内部節点の対角成分には2要素分の寄与が重なります。

コード例2: 全体剛性行列のアセンブリ

import numpy as np

def assemble_global_stiffness(n_elem):
    """一様メッシュ [0,1] の全体剛性行列を組み立てる"""
    n_nodes = n_elem + 1
    h = 1.0 / n_elem
    K = np.zeros((n_nodes, n_nodes))
    ke = (1.0 / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])
    for e in range(n_elem):        # 要素ループ
        for a in range(2):
            for b in range(2):
                K[e + a, e + b] += ke[a, b]   # 対応位置へ加算
    return K

K = assemble_global_stiffness(4)
print("=== 全体剛性行列 (n_elem=4, h=0.25, 境界処理前) ===")
print(K)

# Dirichlet 境界条件: 端の自由度 0 と n を除去する
n_nodes = K.shape[0]
interior = np.arange(1, n_nodes - 1)
print()
print("=== 内部自由度のみの縮約系 (境界処理後) ===")
print(K[np.ix_(interior, interior)])

=== 全体剛性行列 (n_elem=4, h=0.25, 境界処理前) === [[ 4. -4. 0. 0. 0.] [-4. 8. -4. 0. 0.] [ 0. -4. 8. -4. 0.] [ 0. 0. -4. 8. -4.] [ 0. 0. 0. -4. 4.]] === 内部自由度のみの縮約系 (境界処理後) === [[ 8. -4. 0.] [-4. 8. -4.] [ 0. -4. 8.]]

内部節点の対角成分が \( 8 = 2/h \) となり、2要素分が重なっていることが読み取れます。行列は3重対角の対称行列で、有限差分法で現れたラプラシアン行列とよく似た構造を持ちます。

2種類の境界条件

Dirichlet境界条件(Dirichlet condition): 境界での値そのものを指定します(本章では \( u(0)=u(1)=0 \))。解の候補と試験関数の両方に強制する本質的境界条件であり、対応する自由度を系から取り除く(あるいは値を代入する)ことで課します。
Neumann境界条件(Neumann condition): 境界での微分(流束) \( u'(x) \) を指定します。これは弱形式の部分積分で現れた境界項 \( [u'v]_0^1 \) を通じて荷重ベクトルに自然に組み込まれるため、自然境界条件と呼ばれます。 \( u'=0 \) の同次Neumann条件では境界項が消え、特別な処理は不要です。

Dirichlet条件で内部自由度だけの縮約系 \( K_{\text{int}}\,u_{\text{int}} = F_{\text{int}} \) を作れば、あとはこれを解くだけです。系は対称正定値なので、 numpy.linalg.solve で安定に求解できます。

3.5 Python実践: 1次元FEMの実装

ここまでの要素剛性行列・アセンブリ・境界条件処理を1つにまとめ、 1次元Poisson方程式を解く有限要素法ソルバーをNumPyだけでゼロから実装します。荷重ベクトル \( F_i = \int f\,\phi_i\,dx \) は要素ごとに 2点Gauss求積(Gauss quadrature)で数値積分します。 2点Gauss則は3次多項式まで厳密に積分できるため、区分1次要素の荷重計算には十分です。

検証には製作解の方法(method of manufactured solutions)を使います。厳密解を \( u(x)=\sin(\pi x) \) と決め打ちすると、 \( -u''(x) = \pi^2 \sin(\pi x) \) なので \( f(x)=\pi^2\sin(\pi x) \) となります。この \( f \) を入力として数値解を求め、既知の厳密解と比較します。

コード例3: 1次元Poisson方程式のFEMソルバーと厳密解との比較

import numpy as np

def solve_poisson_fem(n_elem, f_func):
    """-u''(x) = f(x) on [0,1], u(0)=u(1)=0 を線形FEMで解く"""
    n_nodes = n_elem + 1
    nodes = np.linspace(0.0, 1.0, n_nodes)
    K = np.zeros((n_nodes, n_nodes))
    F = np.zeros(n_nodes)

    # 基準要素 [-1, 1] 上の2点Gauss求積
    gauss_pts = np.array([-1.0 / np.sqrt(3.0), 1.0 / np.sqrt(3.0)])
    gauss_wts = np.array([1.0, 1.0])

    for e in range(n_elem):
        x_left, x_right = nodes[e], nodes[e + 1]
        h = x_right - x_left
        # 要素剛性行列
        ke = (1.0 / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])
        # 要素荷重ベクトル (Gauss求積)
        fe = np.zeros(2)
        for xi, w in zip(gauss_pts, gauss_wts):
            N = np.array([(1.0 - xi) / 2.0, (1.0 + xi) / 2.0])
            x_phys = x_left + (xi + 1.0) / 2.0 * h   # 物理座標へ写像
            fe += w * (h / 2.0) * f_func(x_phys) * N
        # アセンブリ
        for a in range(2):
            F[e + a] += fe[a]
            for b in range(2):
                K[e + a, e + b] += ke[a, b]

    # Dirichlet 境界条件: 端の自由度を除去
    interior = np.arange(1, n_nodes - 1)
    K_int = K[np.ix_(interior, interior)]
    F_int = F[interior]
    u = np.zeros(n_nodes)
    u[interior] = np.linalg.solve(K_int, F_int)   # 縮約系を求解
    return nodes, u

# 製作解: u(x) = sin(pi x) -> f(x) = pi^2 sin(pi x)
f = lambda x: np.pi**2 * np.sin(np.pi * x)
u_exact = lambda x: np.sin(np.pi * x)

nodes, u = solve_poisson_fem(8, f)
print("nodes :", np.round(nodes, 4))
print("u_FEM :", np.round(u, 6))
print("u_exact:", np.round(u_exact(nodes), 6))

rms = np.sqrt(np.mean((u - u_exact(nodes))**2))
mx = np.max(np.abs(u - u_exact(nodes)))
print("RMS nodal error (n_elem=8): {:.3e}".format(rms))
print("max nodal error (n_elem=8): {:.3e}".format(mx))

nodes : [0. 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1. ] u_FEM : [0. 0.38269 0.707119 0.923895 1.000017 0.923895 0.707119 0.38269 0. ] u_exact: [0. 0.382683 0.707107 0.92388 1. 0.92388 0.707107 0.382683 0. ] RMS nodal error (n_elem=8): 1.110e-05 max nodal error (n_elem=8): 1.665e-05

わずか8要素でも、節点値は厳密解と小数第4位まで一致しています。節点での誤差が \( 10^{-5} \) 程度と極端に小さいのは、 1次元Poisson問題では線形要素の有限要素解が節点上で厳密解に一致する 超収束(superconvergence)という性質によるものです。節点間の内部では誤差はこれより大きくなります。

要素数を増やしたときに、要素内も含めた全体の誤差がどう減るかを調べます。厳密解との差を \( L^2 \) ノルム \( \|u_h - u\|_{L^2} = \left(\int_0^1 (u_h-u)^2\,dx\right)^{1/2} \) で測り、要素幅 \( h \) を半分にするたびに誤差が何倍になるかから収束次数を推定します。

コード例4: 収束性の検証(要素数を変えた誤差評価)

import numpy as np

def l2_error(nodes, u, u_exact_func):
    """各要素で3点Gauss求積により L2 誤差を計算する"""
    gp = np.array([-np.sqrt(3.0/5.0), 0.0, np.sqrt(3.0/5.0)])
    gw = np.array([5.0/9.0, 8.0/9.0, 5.0/9.0])
    err2 = 0.0
    for e in range(len(nodes) - 1):
        xl, xr = nodes[e], nodes[e + 1]
        h = xr - xl
        for xi, w in zip(gp, gw):
            N = np.array([(1.0 - xi)/2.0, (1.0 + xi)/2.0])
            x_phys = xl + (xi + 1.0)/2.0 * h
            uh = N[0]*u[e] + N[1]*u[e + 1]   # 要素内の補間値
            err2 += w * (h/2.0) * (uh - u_exact_func(x_phys))**2
    return np.sqrt(err2)

print("{:>8} {:>10} {:>14} {:>8}".format("n_elem", "h", "L2_error", "rate"))
prev_err, prev_h = None, None
for n_elem in [4, 8, 16, 32, 64]:
    nodes, u = solve_poisson_fem(n_elem, f)
    h = 1.0 / n_elem
    e = l2_error(nodes, u, u_exact)
    if prev_err is None:
        rate = "  -  "
    else:
        rate = "{:.2f}".format(np.log(e/prev_err) / np.log(h/prev_h))
    print("{:>8d} {:>10.5f} {:>14.4e} {:>8}".format(n_elem, h, e, rate))
    prev_err, prev_h = e, h

n_elem h L2_error rate 4 0.25000 3.9127e-02 - 8 0.12500 9.9108e-03 1.98 16 0.06250 2.4859e-03 2.00 32 0.03125 6.2198e-04 2.00 64 0.01562 1.5553e-04 2.00

要素幅 \( h \) を半分にするたびに \( L^2 \) 誤差がおよそ \( 1/4 \) になり、収束次数が2に収束しています。これは区分1次要素の理論値 \( \|u_h - u\|_{L^2} = O(h^2) \) と一致しており、実装が正しいことの強い裏付けになります。

実装のポイント: 要素ループの中で「局所行列を作る → 全体へ加算する」というパターンは、2次元・3次元や高次要素へ拡張しても変わりません。有限要素法の実装は、この局所から全体への組み立ての繰り返しとして理解すると見通しが良くなります。

演習問題

演習3.1: 弱形式の導出

微分方程式 \( -u''(x) + u(x) = f(x) \)、\( x \in (0,1) \)、\( u(0)=u(1)=0 \) について、試験関数 \( v \) を用いた弱形式を部分積分により導出しなさい。本章の \( -u''=f \) と比べて、剛性行列にどのような項が追加されるかを説明しなさい。

演習3.2: Neumann境界条件の実装

コード例3のソルバーを改造し、右端の境界条件を \( u'(1)=g \)(Neumann条件)に変更しなさい。弱形式の境界項 \( [u'v]_0^1 \) が荷重ベクトルにどう寄与するかを示し、 \( g=0 \) のときに端の自由度を除去しなくてよい理由を述べなさい。

演習3.3: 別の製作解での検証

厳密解を \( u(x) = x(1-x) \) とすると、対応する \( f(x) \) はいくつになりますか。その \( f \) をコード例3・4に入力して数値解を求め、\( L^2 \) 誤差を評価しなさい。誤差が要素数によらず機械精度程度に小さくなる理由を、線形要素と2次多項式の関係から考察しなさい(ヒント: 荷重の求積次数にも注意)。

学習目標の確認

本章の学習目標に対する到達度を確認しましょう。

有限要素法が弱形式に基づき、要素と形状関数で解を表現する点を、有限差分法と対比して説明できる
1次元Poisson方程式の弱形式を、部分積分と境界項の消去によって導出できる
線形形状関数から要素剛性行列 \( \frac{1}{h}\begin{pmatrix}1 & -1 \\ -1 & 1\end{pmatrix} \) を導き、全体行列へアセンブリできる
Dirichlet条件(本質的)とNeumann条件(自然)の扱いの違いを説明できる
NumPyだけでFEMソルバーを実装し、収束次数 \( O(h^2) \) を数値的に確認できる

まとめ

有限要素法は、微分方程式を弱形式に緩め、領域を要素に分割して区分的な形状関数で解を近似する手法です。
弱形式は部分積分から得られ、必要な滑らかさが1階微分まで下がるため、折れ線状の近似が可能になります。
区分1次形状関数からは、幅 \( h \) に対して要素剛性行列 \( \frac{1}{h}\begin{pmatrix}1 & -1 \\ -1 & 1\end{pmatrix} \) が導かれます。
要素剛性行列を節点位置で足し合わせるアセンブリにより、3重対角の全体剛性行列が構成されます。
Dirichlet条件は自由度の除去で強制し、Neumann条件は荷重ベクトルに自然に組み込まれます。
製作解の方法で検証すると、線形要素の \( L^2 \) 誤差は理論どおり \( O(h^2) \) で収束します。

次のステップ

本章では、格子点上の差分に依らず、積分(弱形式)を通じて偏微分方程式を解く有限要素法を学びました。第4章「スペクトル法とモンテカルロ法」では、さらに視点を変えて、解を大域的な基底関数(三角関数など)で展開するスペクトル法と、乱数を用いて偏微分方程式や積分を推定するモンテカルロ法を扱います。有限差分法(第1〜2章)・有限要素法(第3章)・スペクトル法(第4章)を並べて眺めることで、数値解法の設計思想の広がりが見えてくるはずです。

参考文献

菊地文雄『有限要素法概説(新訂版)』サイエンス社、1999年.
O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann, 2013.
C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Dover, 2009.
登坂宣好・大西和榮『偏微分方程式の数値解法(第2版)』東京大学出版会、2003年.