有限要素法入門一次元振動問題もっとも単純な問題

有限要素法（Finite Element Method: FEM）は，連続体の力学，熱伝導，電磁気学などの問題を数値的に解くための手法です．基本的な考え方は，複雑な領域を小さな要素に分割し，各要素内で簡単な関数を用いて近似することで，微分方程式を連立一次方程式，非線形方程式，あるいは固有値問題などの代数的な問題に変換することです．

有限要素法の利点は

• 複雑な形状や境界条件に対応しやすい
• 局所的な性質（材料特性など）を容易に取り入れられる
• 汎用性が高く，様々な物理現象に適用可能

です．このページでは，もっとも単純な一次元問題を通じて，有限要素法の基本を学びます．

ここでは，図1 に示す長さ \(L\) の弦の振動問題を考えます．これは，両端が固定された弦に生じる定在波の問題です．この定在波は微分方程式:

\begin{align} \ddiffB{2}{u}{x} + k^2 u = 0 \quad (0 < x < L) \label{eq:vibration_eq} \end{align}

で表されます．$u(x)$ は固有関数，$k$ は波数です．固有値問題としては $λ=k^2$ を固有値とみなします．境界条件は

\begin{align} u(0) = 0, \quad u(L) = 0 \label{eq:bc_vib} \end{align}

です．

ここで，問題を整理します．解くべきことは，

です．これは，微分方程式の固有値問題の一例です．ここで用いた例は非常に単純なため，解析解があります．複雑な問題，例えば境界が複雑な二次元や三次元問題になると，解析解がありません．そうなると，有限要素法を使って数値的に解くことが行われます．

ここで設定した問題の解析解は，\(u_n(x) = \sin\left(n\pi x/L\right)\)， \(k_n = n\pi/L\) です（\(n = 1, 2, 3, \dots\)）．解析解があるので，有限要素法の結果と比較して，数値解の精度を評価できます．

 注意 (物理) 振動の問題というからには，振動数が欲しいところです．この弦の振動の波の速度を $c$ とすると，$c=\omega_n/k_n$ なので，$\omega_n = c k_n$ です $(n = 1,\,2,\,3,\,\dots)$．力学から，速度は $c=\sqrt{T/\rho}$ で，$T$ は張力，$\rho$ は線密度です．要するに波数: $k_n$ と弦に関わる特性 (張力: $T$，線密度: $\rho$) が分かれば，振動数が分かるということです．

有限要素法では，微分方程式を弱形式（variational form）に変換します．具体的には，式\eqref{eq:vibration_eq} の両辺にテスト関数 \(v(x)\) を掛け，区間 \([0,L]\) で積分します．ここで，テスト関数 \(v(x)\) は次の条件を満たすものとします．

• 区間 \([0, L]\) で十分滑らか，
• 両端固定条件に対応して \(v(0) = v(L) = 0\) を満たす，

すると，

\begin{align} \int_0^L \ddiffB{2}{u}{x} v \, \diff x + k^2 \int_0^L u v \, \diff x = 0 \label{eq:weak1_vib} \end{align}

が得られます．テスト関数 \(v(x)\) が区間 \([0, L]\) で正則であれば，弱形式の式\eqref{eq:weak1_vib}は微分方程式\eqref{eq:vibration_eq}を満たすための必要条件です．弱形式の式\eqref{eq:weak1_vib}がいかなる正則なテスト関数: $v(x)$で成立するならば，弱形式の式\eqref{eq:weak1_vib}は微分方程式\eqref{eq:vibration_eq}を満たす必要十分条件になります．

テスト関数を使う弱形式の利点は以下の通りです．

• Neumann 型境界条件を自然境界条件として扱いやすい．
• 古典解（強解）より広い解のクラス（弱解）を許容し，数値計算に適した形にする．
• Galerkin法 (後で説明) により，近似関数との内積を取って離散化を行う基盤となる．

式\eqref{eq:weak1_vib}の左辺第一項を部分積分すると，

\begin{align} \left[ \ddiffA{u}{x} v \right]_0^L - \int_0^L \ddiffA{u}{x} \ddiffA{v}{x} \, \diff x + k^2 \int_0^L u v \, \diff x = 0 \label{eq:weak2_vib} \end{align}

となります．両端固定の境界条件は本質境界条件として近似解に課します．そのため，テスト関数は $v(0)=v(L)=0$ を満たす関数から選びます．したがって境界項 $[u^\prime v]_0^L$ は 0 になります．すると，

\begin{align} \int_0^L \ddiffA{u}{x} \ddiffA{v}{x} \, \diff x = k^2 \int_0^L u v \, \diff x \label{eq:weak_form_vib} \end{align}

です．これが弱形式です．固有値問題として，固有値: \(k^2\) と固有関数: \(u(x)\) を求めることになります．

 Galerkin 法  有限要素法では，弱形式をそのまま計算するのではなく，有限個の関数で近似します．このとき，未知関数 \(u(x)\) を表すための近似関数と，弱形式に現れるテスト関数 \(v(x)\) を同じ形状関数から選ぶ方法を Galerkin 法と呼びます．たとえば \(u(x)\) を節点値と形状関数の線形結合で表すなら，テスト関数にも各節点に対応する形状関数を用います．

Galerkin 法では，近似解が微分方程式を点ごとに完全に満たすのではなく，残差が選ばれたテスト関数すべてに対して直交するように求められます．これにより，連続的な微分方程式を有限個の未知数に対する代数方程式へ変換できます．

有限要素法では，弱形式の式\eqref{eq:weak_form_vib}を計算します．式中の \(u(x)\) や \(v(x)\) は連続量なので，数値計算では扱いが難しいです．そこで，コンピューターでの計算を容易にするために，これらを離散化します．そのために，領域 \([0, L]\) を \(n\) 個の要素に分割し，各要素 \(e\) 内で解 \(u\) を線形関数で近似:

\begin{align} u^e(x) = a^e + b^e x \label{eq:approx} \end{align}

します (図2)．ここで，上付きの \(e\) は「要素番号」や「要素ごとの局所量」を示します．この式は，まだ扱いにくいので，節点値と形状関数を用いて表現します．具体的には，

\begin{align} u^e(x) = u_i^e \phi_i^e(x) + u_{i+1}^e \psi_i^e(x) \label{eq:shape_approx} \end{align}

です (図3)．ここで，\(u^e(x)\) は要素 \(e\) 内で定義された近似解，\(u_i^e\) は要素 \(e\) の節点 \(i\) における節点値です．要素の定義域は \([x_i, x_{i+1}]\) です．ここで，$(\phi_i^e,\,\psi_i^e )$ は要素 \(e\) に対応する形状関数です．具体的には，

\begin{align} \phi_i^e(x) = \frac{x_{i+1} - x}{h^e}, \quad \psi_i^e(x) = \frac{x - x_i}{h^e} \label{eq:shape_functions} \end{align}

です．ここで，\(h^e = x_{i+1} - x_i\) は要素長です．これらの形状関数の特徴は

• \(\phi_i^e(x)\)：節点 \(i\) で1，節点 \(i+1\) で0となり，要素内で線形に変化する．
• \(\psi_i^e(x)\)：節点 \(i\) で0，節点 \(i+1\) で1となり，要素内で線形に変化する．

です．これらの形状関数を用いることで，要素内での近似関数は節点値と形状関数の組み合わせで簡潔に表現できます．

領域 $[0, L]$ での弱形式 \eqref{eq:weak_form_vib} を，領域 $[x_i, x_{i+1}] $ の要素で弱形式を適用します．するとこの要素 \(e\) での弱形式は

\begin{align} \int_{x_i}^{x_{i+1}} \ddiffA{u^e}{x} \ddiffA{v^e}{x} \, \diff x = k^2 \int_{x_i}^{x_{i+1}} u^e v^e \, \diff x \label{eq:element_weak_vib} \end{align}

です．これを，Galerkin 法を用いて解きます．すなわち，テスト関数は形状関数: $v^e=\phi_i^e(x)$ と $v^e=\psi_i^e(x)$ を使います．

 テスト関数: $\phi_i^e(x)$ の場合  式\eqref{eq:element_weak_vib} にテスト関数: $v=\phi_i^e(x)$ を使います．解の関数は式\eqref{eq:shape_approx}です．これらから，

\begin{align} \ddiffA{u^e}{x} = u_i^e \ddiffA{\phi_i^e}{x} + u_{i+1}^e \ddiffA{\psi_i^e}{x} , \quad \ddiffA{v^e}{x}=\ddiffA{\phi_i^e}{x} \label{eq:element_weak_phi_vib} \end{align}

が得られます．これらを式\eqref{eq:element_weak_vib}に代入して積分すると，

\begin{align} \int_{x_i}^{x_{i+1}} \left( \ddiffA{\phi_i^e}{x}\ddiffA{\phi_i^e}{x}u_i^e + \ddiffA{\phi_i^e}{x}\ddiffA{\psi_i^e}{x}u_{i+1}^e \right)\,\diff x = k^2 \int_{x_i}^{x_{i+1}} \left(\phi_i^e\phi_i^eu_i^e + \phi_i^e\psi_i^e u_{i+1}^e \right) \, \diff x \label{eq:element_phi_vib} \end{align}

となります．

 テスト関数: $\psi_i^e(x)$ の場合  同様に，式\eqref{eq:element_weak_vib} にテスト関数: $v=\psi_i^e(x)$ を使います．解の関数は式\eqref{eq:shape_approx}です．先と同じ計算を行うと，

\begin{align} \int_{x_i}^{x_{i+1}} \left( \ddiffA{\psi_i^e}{x}\ddiffA{\phi_i^e}{x}u_i^e + \ddiffA{\psi_i^e}{x}\ddiffA{\psi_i^e}{x}u_{i+1}^e \right)\,\diff x = k^2 \int_{x_i}^{x_{i+1}} \left(\psi_i^e\phi_i^eu_i^e + \psi_i^e\psi_i^e u_{i+1}^e \right) \,\diff x \label{eq:element_psi_vib} \end{align}

となります．

 行列形式に変換  式\eqref{eq:element_phi_vib} と式\eqref{eq:element_psi_vib} をまとめると，

\begin{align} \begin{bmatrix} \int_{x_i}^{x_{i+1}}\ddiffA{\phi_i^e}{x}\ddiffA{\phi_i^e}{x}\diff x & \int_{x_i}^{x_{i+1}}\ddiffA{\phi_i^e}{x}\ddiffA{\psi_i^e}{x}\diff x \\ \int_{x_i}^{x_{i+1}}\ddiffA{\psi_i^e}{x}\ddiffA{\phi_i^e}{x}\diff x & \int_{x_i}^{x_{i+1}}\ddiffA{\psi_i^e}{x}\ddiffA{\psi_i^e}{x}\diff x \end{bmatrix} \begin{bmatrix} u_i^e \\ u_{i+1}^e \end{bmatrix} = k^2 \begin{bmatrix} \int_{x_i}^{x_{i+1}} \phi_i^e\phi_i^e\diff x & \int_{x_i}^{x_{i+1}} \phi_i^e\psi_i^e\diff x \\ \int_{x_i}^{x_{i+1}} \psi_i^e\phi_i^e\diff x & \int_{x_i}^{x_{i+1}} \psi_i^e\psi_i^e\diff x \end{bmatrix} \begin{bmatrix} u_i^e \\ u_{i+1}^e \end{bmatrix} \label{eq:element_matrix_vib} \end{align}

となります．それぞれの積分を計算すると，

\begin{align} \int_{x_i}^{x_{i+1}}\ddiffA{\phi_i^e}{x}\ddiffA{\phi_i^e}{x}\diff x = \frac{1}{h^e}, \quad \int_{x_i}^{x_{i+1}}\ddiffA{\phi_i^e}{x}\ddiffA{\psi_i^e}{x}\diff x = -\frac{1}{h^e}, \quad \int_{x_i}^{x_{i+1}}\ddiffA{\psi_i^e}{x}\ddiffA{\phi_i^e}{x}\diff x = -\frac{1}{h^e}, \quad \int_{x_i}^{x_{i+1}}\ddiffA{\psi_i^e}{x}\ddiffA{\psi_i^e}{x}\diff x = \frac{1}{h^e}\\[10pt] \int_{x_i}^{x_{i+1}} \phi_i^e\phi_i^e\diff x = \frac{h^e}{3}, \quad \int_{x_i}^{x_{i+1}} \phi_i^e\psi_i^e\diff x = \frac{h^e}{6}, \quad \int_{x_i}^{x_{i+1}} \psi_i^e\phi_i^e\diff x = \frac{h^e}{6}, \quad \int_{x_i}^{x_{i+1}} \psi_i^e\psi_i^e\diff x = \frac{h^e}{3} \label{eq:mass_integral_vib} \end{align}

です．この結果を使うと，式\eqref{eq:element_matrix_vib}を行列形式で表すと，

\begin{align} \mathbf{k}^e \mathbf{u}^e = k^2 \mathbf{m}^e \mathbf{u}^e \label{eq:element_eq_vib} \end{align}

となります．ここで，要素剛性マトリックス \(\mathbf{k}^e\) と要素質量マトリックス \(\mathbf{m}^e\) は，以下の通りです．

\begin{align} \mathbf{k}^e = \frac{1}{h^e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \qquad \mathbf{m}^e = \frac{h^e}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \label{eq:elem:km_matrix} \end{align}

前節では，一つの要素 $e$ に対して，式 \eqref{eq:element_weak_vib} から出発し，要素方程式である式\eqref{eq:element_eq_vib}を得ました．要素方程式の行列は式 \eqref{eq:elem:km_matrix} で与えられます．ここでは，この結果を使い，全体を表す方程式を組み立てます．

区間 \([0,\ L]\) を \(N-1\) 個の要素に分割します．その様子を図4に示します．この図から，

• 節点の座標は，$x_1,\ x_2,\ x_3,\ \dots,\ x_N$ です．
• 各要素の間隔は，$h_1,\ h_2,\ h_3,\ \dots,\ h_{N-1}$ です．
• 領域全体の節点値ベクトルは \(\mathbf{u} = [u_1,\ u_2,\ u_3,\ \dots,\ u_{N}]^T\) です．

です．この分割のもと，式 \eqref{eq:weak_form_vib} を計算します．この式は，

\begin{align} \int_{x_1}^{x_{2}} \ddiffA{u_1^e}{x} \ddiffA{v_1^e}{x} \, \diff x + &\int_{x_2}^{x_{3}} \ddiffA{u_2^e}{x} \ddiffA{v_2^e}{x} \, \diff x + \int_{x_3}^{x_{4}} \ddiffA{u_3^e}{x} \ddiffA{v_3^e}{x} \, \diff x + \cdots + \int_{x_{N-1}}^{x_N} \ddiffA{u_{N-1}^e}{x} \ddiffA{v_{N-1}^e}{x} \, \diff x = \\[10pt] &k^2 \int_{x_1}^{x_2} u_1^e v_1^e \, \diff x + k^2 \int_{x_2}^{x_3} u_2^e v_2^e \, \diff x + k^2 \int_{x_3}^{x_4} u_3^e v_3^e \, \diff x + \cdots + k^2 \int_{x_{N-1}}^{x_N} u_{N-1}^e v_{N-1}^e \, \diff x \end{align}

となります．ここで，各要素内の近似関数: $u^e$ とテスト関数: $v^e$ には，前節と同じ形状関数を使います．すると，この式の右辺と左辺のそれぞれの項は，式\eqref{eq:element_matrix_vib} から式\eqref{eq:elem:km_matrix} を用いて計算できます．全体剛性マトリックスは

\begin{align} \mathbf{K} = \begin{bmatrix} 1/h_1 & -1/h_1 & 0 & \dots & 0 \\ -1/h_1 & 1/h_1+1/h_2 & -1/h_2 & \dots & 0\\ 0 & -1/h_2 & 1/h_2+1/h_3 & \dots & \vdots \\ \vdots & \vdots & \vdots & \ddots & -1/h_{N-1}\\ 0 & 0 & \vdots & -1/h_{N-1} & 1/h_{N-1} \end{bmatrix} \end{align}

です．各成分は，

\begin{align} &K_{i,i} = \frac{1}{h_{i-1}} + \frac{1}{h_i},& &K_{i,i+1} = K_{i+1,i} = -\frac{1}{h_i} \quad (i=2,3,\dots,N-2)&\\[10pt] &K_{1,1} = \frac{1}{h_1},& &K_{N,N} = \frac{1}{h_{N-1}}& \end{align}

です．

同様に，全体質量マトリックスは

\begin{align} \mathbf{M} = \begin{bmatrix} h_1/3 & h_1/6 & 0 & \dots & 0 \\ h_1/6 & (h_1+h_2)/3 & h_2/6 & \dots & 0 \\ 0 & h_2/6 & (h_2+h_3)/3 & \dots & \vdots \\ \vdots & \vdots & \vdots & \ddots & h_{N-1}/6\\ 0 & 0 & \dots & h_{N-1}/6 & h_{N-1}/3 \end{bmatrix} \end{align}

となります．各成分は，

\begin{align} &M_{i,i} = \frac{(h_i+h_{i+1})}{3},& &M_{i,i+1} = M_{i+1,i} = -\frac{h_i}{6} \quad (i=2,3,\dots,N-2)&\\[10pt] &M_{1,1} = \frac{h_1}{3},& &M_{N,N} = \frac{h_{N-1}}{3}& \end{align}

です．

これらのマトリックスは，いずれも $N\times N$ の三重対角の対称行列 (エルミート行列) です．

以上により，領域全体での離散化された固有値問題

\begin{align} \mathbf{K}\mathbf{u}=k^2 \mathbf{M}\mathbf{u} \end{align}

が得られます．ここで，

\begin{align} \mathbf{u}= \begin{bmatrix} u_1 & u_2 & \cdots & u_N \end{bmatrix}^{\mathrm{T}} \end{align}

です．

 境界条件  両端固定の境界条件は

\begin{align} u_1 = 0, \qquad u_N = 0 \end{align}

です．したがって，未知数として残るのは内部節点の値

\begin{align} u_2,\ u_3,\ \dots,\ u_{N-1} \end{align}

です．境界条件を適用するには，全体行列 \(\mathbf{K}\)，\(\mathbf{M}\) から 1 行目・1 列目と \(N\) 行目・\(N\) 列目を取り除き，内部節点だけの行列を作ります．これにより，

\begin{align} \mathbf{K}_{\mathrm{in}}\mathbf{u}_{\mathrm{in}}= k^2 \mathbf{M}_{\mathrm{in}}\mathbf{u}_{\mathrm{in}} \end{align}

という一般化固有値問題が得られます．この固有値問題を解くことで，固有値 \(k^2\) と，内部節点での固有モード \(\mathbf{u}_{\mathrm{in}}\) が求まります．境界節点の値 \(u_1=0\)，\(u_N=0\) を加えれば，区間全体での振動モードが得られます．

 固有ベクトル・固有値の性質  $\mathbf{K}$ と $\mathbf{M}$ は，それぞれ対称行列です．境界条件を適用した後のは $\mathbf{K}_{\mathrm{in}}$ と $\mathbf{M}_{\mathrm{in}}$ であり，これらは正定値です．したがって，一般化固有値問題の固有値: $k^2$ は実数で，正の値になります．固有ベクトル: $\mathbf{u}_{\mathrm{in}}$ も正規直交になります．

プログラム例

この節では，前節までの理論を実装した数値例を示します．以下のコード (fem1d_vibration.py) は，Python を使って有限要素法の固有値問題を解く例です．実行方法は，

> python  fem1d_vibration.py 100  3

です．ここで，コマンドライン引数の 100 は要素数を表し，3 は求める固有モードの数を表します．このプログラムは，両端固定の弦の振動モードを有限要素法で計算し，最初の 3 つの固有モードをファイル (fem1d_results.txt) に保存します．

            #!/usr/bin/env python3
"""Finite element solution of a 1D string vibration eigenvalue problem.

Problem:
    u'' + k^2 u = 0,  0 < x < L
    u(0) = u(L) = 0

Weak form:
    integral_0^L u' v' dx = k^2 integral_0^L u v dx

This script uses linear finite elements on a uniform mesh.  The command-line
argument N is the number of elements, so the number of nodes is N + 1.
The boundary nodes are fixed, and the generalized eigenvalue problem is solved
for the internal degrees of freedom.
"""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np


LENGTH = 1.0
OUTPUT_FILE = "fem1d_results.txt"


def build_global_matrices(num_elements: int, length: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Build global stiffness and mass matrices for linear 1D elements."""
    num_nodes = num_elements + 1
    nodes = np.linspace(0.0, length, num_nodes)

    stiffness = np.zeros((num_nodes, num_nodes))
    mass = np.zeros((num_nodes, num_nodes))

    for elem in range(num_elements):
        h = nodes[elem + 1] - nodes[elem]

        elem_stiffness = (1.0 / h) * np.array(
            [
                [1.0, -1.0],
                [-1.0, 1.0],
            ]
        )
        elem_mass = (h / 6.0) * np.array(
            [
                [2.0, 1.0],
                [1.0, 2.0],
            ]
        )

        global_nodes = [elem, elem + 1]
        for local_i, global_i in enumerate(global_nodes):
            for local_j, global_j in enumerate(global_nodes):
                stiffness[global_i, global_j] += elem_stiffness[local_i, local_j]
                mass[global_i, global_j] += elem_mass[local_i, local_j]

    return nodes, stiffness, mass


def apply_fixed_end_boundary_conditions(
    stiffness: np.ndarray,
    mass: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    """Remove both fixed boundary degrees of freedom."""
    return stiffness[1:-1, 1:-1], mass[1:-1, 1:-1]


def solve_generalized_symmetric_eigenproblem(
    stiffness: np.ndarray,
    mass: np.ndarray,
    num_modes: int,
) -> tuple[np.ndarray, np.ndarray]:
    """Solve K u = lambda M u using only NumPy.

    Since M is symmetric positive definite, write M = L L^T and transform
    the problem into a standard symmetric eigenvalue problem:

        L^{-1} K L^{-T} y = lambda y,  u = L^{-T} y.
    """

    cholesky = np.linalg.cholesky(mass)
    temp = np.linalg.solve(cholesky, stiffness)
    standard_matrix = np.linalg.solve(cholesky, temp.T).T

    eigenvalues, transformed_vectors = np.linalg.eigh(standard_matrix)
    order = np.argsort(eigenvalues)
    eigenvalues = eigenvalues[order][:num_modes]
    transformed_vectors = transformed_vectors[:, order][:, :num_modes]

    eigenvectors = np.linalg.solve(cholesky.T, transformed_vectors)

    for mode in range(num_modes):
        norm = np.sqrt(eigenvectors[:, mode].T @ mass @ eigenvectors[:, mode])
        eigenvectors[:, mode] /= norm

    return eigenvalues, eigenvectors


def expand_internal_modes(eigenvectors_internal: np.ndarray, num_elements: int) -> np.ndarray:
    """Add zero boundary values to internal eigenvectors."""
    num_nodes = num_elements + 1
    full_vectors = np.zeros((num_nodes, eigenvectors_internal.shape[1]))
    full_vectors[1:-1, :] = eigenvectors_internal
    return full_vectors


def write_results(
    path: Path,
    nodes: np.ndarray,
    eigenvalues: np.ndarray,
    eigenvectors: np.ndarray,
) -> None:
    """Write eigenvalues and eigenvectors to a text file."""
    with path.open("w", encoding="utf-8") as file:
        file.write("# 1D FEM string vibration eigenvalue problem\n")
        file.write("# Equation: u'' + k^2 u = 0, u(0)=u(L)=0\n")
        file.write("# Eigenvalue lambda = k^2\n")
        file.write(f"# number_of_elements = {len(nodes) - 1}\n")
        file.write(f"# number_of_nodes = {len(nodes)}\n")
        file.write(f"# length = {nodes[-1] - nodes[0]:.16e}\n\n")

        file.write("# Eigenvalues\n")
        file.write("# mode  lambda=k^2              k\n")
        for index, eigenvalue in enumerate(eigenvalues, start=1):
            file.write(f"{index:5d}  {eigenvalue:.16e}  {np.sqrt(eigenvalue):.16e}\n")

        file.write("\n# Eigenvectors\n")
        file.write("# Each eigenvector includes the fixed boundary values u(0)=u(L)=0.\n")
        header = "node_index  x" + "".join(f"  mode_{i}" for i in range(1, len(eigenvalues) + 1))
        file.write(f"# {header}\n")
        for node_index, x in enumerate(nodes, start=1):
            values = "".join(f"  {eigenvectors[node_index - 1, mode]: .16e}" for mode in range(len(eigenvalues)))
            file.write(f"{node_index:10d}  {x:.16e}{values}\n")


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Solve the 1D FEM string vibration eigenvalue problem."
    )
    parser.add_argument("N", type=int, help="number of elements")
    parser.add_argument("m", type=int, help="number of modes, counted from the smallest eigenvalue")
    return parser.parse_args()


def main() -> None:
    args = parse_args()

    if args.N < 2:
        raise SystemExit("N must be at least 2 because both ends are fixed.")

    max_modes = args.N - 1
    if args.m < 1 or args.m > max_modes:
        raise SystemExit(f"m must satisfy 1 <= m <= {max_modes}.")

    nodes, stiffness, mass = build_global_matrices(args.N, LENGTH)
    stiffness_internal, mass_internal = apply_fixed_end_boundary_conditions(stiffness, mass)
    eigenvalues, eigenvectors_internal = solve_generalized_symmetric_eigenproblem(
        stiffness_internal,
        mass_internal,
        args.m,
    )
    eigenvectors = expand_internal_modes(eigenvectors_internal, args.N)

    output_path = Path(OUTPUT_FILE)
    write_results(output_path, nodes, eigenvalues, eigenvectors)
    print(f"Wrote {args.m} modes to {output_path}")


if __name__ == "__main__":
    main()
            

一般化固有値問題の解法

 一般化固有値問題の解法  有限要素法により，最終的に次の一般化固有値問題が得られます．

\begin{align} \mathbf{K}\mathbf{u} = \lambda \mathbf{M}\mathbf{u} \label{eq:generalized_eigenproblem} \end{align}

ここで，\(\lambda=k^2\) です．\(\mathbf{K}\) は剛性マトリックス，\(\mathbf{M}\) は質量マトリックスです．両端固定の境界条件を適用した後，\(\mathbf{K}\) と \(\mathbf{M}\) は対称行列になります．また，\(\mathbf{M}\) は正定値行列です．そこで，\(\mathbf{M}\) を Cholesky 分解します．

\begin{align} \mathbf{M}=\mathbf{L}\mathbf{L}^{\mathrm{T}} \label{eq:cholesky_mass} \end{align}

ここで，\(\mathbf{L}\) は下三角行列です．式\eqref{eq:generalized_eigenproblem} に \(\mathbf{M}=\mathbf{L}\mathbf{L}^{\mathrm{T}}\) を代入すると，

\begin{align} \mathbf{K}\mathbf{u}=\lambda \mathbf{L}\mathbf{L}^{\mathrm{T}}\mathbf{u} \end{align}

となります．ここで，新しい未知ベクトル \(\mathbf{y}\) を

\begin{align} \mathbf{y}=\mathbf{L}^{\mathrm{T}}\mathbf{u} \label{eq:y_transform} \end{align}

と定義します．この式から，

\begin{align} \mathbf{u}=\mathbf{L}^{-\mathrm{T}}\mathbf{y} \label{eq:u_transform} \end{align}

です．これを一般化固有値問題に代入すると，

\begin{align} \mathbf{K}\mathbf{L}^{-\mathrm{T}}\mathbf{y}=\lambda \mathbf{L}\mathbf{y} \end{align}

となります．さらに左から \(\mathbf{L}^{-1}\) を掛けると，

\begin{align} \mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\mathrm{T}}\mathbf{y}=\lambda \mathbf{y} \label{eq:standard_eigenproblem} \end{align}

が得られます．これは通常の固有値問題です．したがって，

\begin{align} \mathbf{A}=\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\mathrm{T}} \end{align}

とおけば，

\begin{align} \mathbf{A}\mathbf{y}=\lambda \mathbf{y} \end{align}

を解けばよいことになります．この固有値問題を解いて得られる固有値 \(\lambda\) は，元の一般化固有値問題の固有値と同じです．一方，得られる固有ベクトルは \(\mathbf{y}\) なので，元の節点値ベクトル \(\mathbf{u}\) に戻すには，式\eqref{eq:u_transform} を用いて

\begin{align} \mathbf{u}=\mathbf{L}^{-\mathrm{T}}\mathbf{y} \end{align}

を計算します．

Python プログラムでは，逆行列 \(\mathbf{L}^{-1}\) や \(\mathbf{L}^{-\mathrm{T}}\) を明示的には作りません．数値計算では，逆行列を直接求めるよりも，連立一次方程式として解く方が安定で効率的だからです．そのため，NumPy の np.linalg.solve を使って，

\begin{align} \mathbf{L}\mathbf{X}=\mathbf{K} \end{align}

を解くことで，

\begin{align} \mathbf{X}=\mathbf{L}^{-1}\mathbf{K} \end{align}

を求めます．続いて，

\begin{align} \mathbf{A}=\mathbf{X}\mathbf{L}^{-\mathrm{T}}=\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\mathrm{T}} \end{align}

を作ります．これがプログラム中の standard_matrix です．

このページでは，一次元振動問題を通じて有限要素法の基本を学びました．有限要素法は，加速器のRF共振器設計などで重要な固有値問題を効率的に解くことができます．

より複雑な問題（2次元，3次元，非線形）にも同様の考え方を拡張できます．

参考資料

更新履歴