固有値問題の離散化：差分・有限要素の固有値の収束と上界・下界

0. この記事の要点

弦や太鼓の膜の振動の音の高さ、量子力学のエネルギーの値は、どれも微分作用素の固有値という数で決まります。コンピュータでは、微分作用素を行列に置き換えて、その行列の固有値で近似します。近似値は真の値より大きくなる方法と小さくなる方法があり、誤差の減り方も固有値の大きさによって変わります。この記事では、格子を細かくしたときに誤差がどちらの向きに、どれだけの速さで減るか、下から何番目までの固有値が信頼できるかを調べます。読み終えると、上からの近似と下からの近似で真の値を挟む考え方が分かります。

• 区間の Dirichlet ラプラシアン（両端で 0 になる条件をつけた $-d^2/dx^2$）の 3 点差分の固有値は $\lambda_k^{\mathrm{FD}} = \frac{4}{h^2}\sin^2\frac{k\pi h}{2}$ で、厳密な固有値 $\lambda_k = k^2\pi^2$ に下から近づき、$0 \lt \lambda_k - \lambda_k^{\mathrm{FD}} \le \lambda_k^2h^2/12$ です。誤差は $\lambda_k^2$ に比例するので、信頼できるのは $kh$ が小さい下位の固有値だけです。
• 適合な Galerkin 近似（$V_h \subset V$。近似に使う関数の空間をもとの空間の部分空間にとること）の Ritz 値（近似の行列の固有値）は、厳密な固有値に上から近づきます（$\lambda_k \le \lambda_k^h$）。部分空間を広げると単調に減ります。
• Ritz 値の誤差は、固有関数のエネルギーノルム（微分の大きさまで含めた「距離」）での近似誤差の 2 乗の次数です。Ritz 射影（近似の空間の中でエネルギーノルムの意味で最も近い関数をとる操作。§2 で定義します）の誤差が $\|w - P_hw\|_a \le \eta\|Aw\|$（$A$ は対象の作用素、$w$ はその定義域の関数）を満たせば $\lambda_k^h - \lambda_k \le 4\eta^2\lambda_k^2$ で、区分 1 次の有限要素法では $\eta = O(h)$ なので $O(\lambda_k^2h^2)$ です。
• 1 次元の等間隔の区分 1 次要素では $\lambda_k^h - \lambda_k = \lambda_k^2h^2/12 + O(h^4)$ で、質量行列を集中化する（質量行列を対角行列に置き換える）と 3 点差分に一致して符号が逆になります。両者の平均は $O(h^4)$ です。
• 単純固有値の Ritz ベクトルの誤差は $O(h)$（エネルギーノルム）で、定数は固有値のすき間の逆数に比例します。
• 全空間の問題を有限の箱で切って（両端で 0 とおいて）解くと、近似固有値は上から単調に近づき、固有関数が $e^{-\kappa|x|}$ で減衰するなら誤差は $O(e^{-2\kappa R})$ です。
• 下に有界な自己共役作用素では、本質スペクトル（有限重複度の孤立した固有値以外のスペクトルの部分）の下端より下に偽の固有値（作用素のスペクトルと無関係に行列にだけ現れる固有値）は現れません。すき間の中の固有値にはこの保証がありません。

この記事で使う事前知識

• Hilbert 空間と正規直交基底: 内積をもち、極限をとっても外に出ない空間と、その中の互いに直交する長さ 1 のベクトルからなる基底。→ Hilbert 空間、正規直交基底
• Dirichlet ラプラシアン: 領域の境界で関数が 0 になる条件をつけたラプラシアン。固有値は小さい方から並べられ、無限大に発散します。→ Dirichlet ラプラシアンと Weyl の法則
• ミニマックス値: 固有値を、部分空間を動かして Rayleigh 商の最大値を最小にしたものとして表した量。→ ミニマックス値
• Rayleigh–Ritz の上界: 部分空間に制限して求めた近似固有値は、真の固有値より小さくならないという事実。→ Rayleigh–Ritz の上界
• Galerkin 近似と Ritz 射影: 有限次元の部分空間の中で、エネルギーノルムの意味で最もよく近似する関数をとる操作。→ Galerkin 近似
• 直交射影: ベクトルを部分空間の上に、誤差が部分空間と直交するように落とす操作。→ 直交射影

1. 動機

量子力学のハミルトニアン、膜や梁の振動、拡散の緩和時間。どれも「微分作用素の固有値を求める」問題で、計算では作用素を行列に置き換えて対角化します。行列の固有値を求める算法（QR 法、Lanczos 法、逆反復）の誤差は丸めの大きさまで小さくできますが（対称固有値問題と特異値分解）、得られるのは行列の固有値です。それが作用素の固有値からどれだけ離れているかは、離散化の問題です。

この章で答える問いは 3 つです。

1. 格子幅 $h$ を半分にすると、固有値の誤差はどれだけ減るか（次数）。
2. 近似固有値は真の値より大きいか小さいか（向き）。
3. $n$ 次の行列の $n$ 個の固有値のうち、どこまでが信頼できるか。

2 が分かれば、上からの近似と下からの近似で真の値を挟めます。固有値の精度保証の章は、その挟み込みを丸め誤差まで含めて証明に変える方法を扱っていました。この章はその手前で、離散化の誤差の次数と向きを、差分法・有限要素法・スペクトル法・領域の切断について調べます。Rayleigh–Ritz の上界、Temple の不等式、Lehmann–Goerisch の下界、Liu の下界（どれも、近似から真の固有値の上や下の限界を与える不等式）は、既存の章の定理を使います（定理 3.1［無限次元の固有値の包含］、定理 4.1［無限次元の固有値の包含］、定理 5.2［無限次元の固有値の包含］、定理 6.3［無限次元の固有値の包含］）。

2. 設定：形式の固有値問題と Galerkin 近似

固有値の精度保証の章の §2 と同じ設定をとります。$\mathcal{H}$ を実 Hilbert 空間（内積 $\langle\cdot,\cdot\rangle$、ノルム $\|\cdot\|$）、$A$ を $\mathcal{H}$ 上の自己共役作用素で、固有値 $0 \lt \lambda_1 \le \lambda_2 \le \cdots \to \infty$（重複度こみ）と、固有ベクトルからなる正規直交基底 $\{u_m\}$（$Au_m = \lambda_mu_m$）をもつとします。$w = \sum_mc_mu_m$ について

$$V = \Bigl\{w : \sum_m\lambda_mc_m^2 \lt \infty\Bigr\},\qquad a(w, w') = \sum_m\lambda_mc_mc_m',\qquad \|w\|_a = a(w, w)^{1/2}$$

を $A$ の形式とその定義域、エネルギーノルムとします（微分作用素なら、$a(w,w)$ は $\int|w'|^2$ のように「微分の大きさの 2 乗」を測る量です）。$D(A) = \{w : \sum_m\lambda_m^2c_m^2 \lt \infty\} \subset V$ で、$w \in D(A)$ について $\|Aw\|^2 = \sum_m\lambda_m^2c_m^2$、$a(w, v) = \langle Aw, v\rangle$（$v \in V$）です。とくに $a(u_m, v) = \lambda_m\langle u_m, v\rangle$ です。$A^{-1}$ は $\mathcal{H}$ 全体で定義された有界作用素（$A^{-1}w = \sum_m(c_m/\lambda_m)u_m$）で、値域は $D(A)$ です。$\lambda_1\|w\|^2 \le \|w\|_a^2$ です。

モデル問題は次の 2 つです。

• 区間: $\mathcal{H} = L^2(0,1)$、$V = H^1_0(0,1)$、$a(u,v) = \int_0^1u'v'\,dx$。$A = -d^2/dx^2$（Dirichlet 条件）、$\lambda_k = k^2\pi^2$、$u_k = \sqrt2\sin k\pi x$（命題 6.2［Sturm–Liouville 問題と Green 関数］）、$D(A) = H^2(0,1)\cap H^1_0(0,1)$、$\|Aw\| = \|w''\|$。
• 有界領域 $\Omega$ の Dirichlet ラプラシアン: $V = H^1_0(\Omega)$、$a(u,v) = \int_\Omega\nabla u\cdot\nabla v\,dx$（Dirichlet ラプラシアンと Weyl の法則）。

次の補題は、固有値を「部分空間を動かして Rayleigh 商の最大値を最小にしたもの」として書き直します。後の上界の証明は、動かす部分空間の範囲が狭いほど最小値が大きくなる、という一点に帰着します。

補題 2.1（形式についてのミニマックス原理）✓

上の設定で、$R(w) = a(w,w)/\|w\|^2$（$w \in V$、$w \ne 0$）を Rayleigh 商とすると、$k \ge 1$ について

$$\lambda_k = \min_{\substack{F \subset V\\ \dim F = k}}\ \max_{w \in F,\ w \ne 0}R(w)$$

が成り立つ。

証明（補題 2.1）

$F_0 = \mathrm{span}\{u_1, \ldots, u_k\}$ とすると、$w = \sum_{m\le k}c_mu_m$ について $a(w,w) = \sum_{m\le k}\lambda_mc_m^2 \le \lambda_k\|w\|^2$ で、$w = u_k$ で等号なので、$\max_{F_0}R = \lambda_k$ です。逆に、$F \subset V$ を任意の $k$ 次元部分空間とします。線形写像 $F \to \mathbb{R}^{k-1}$、$w \mapsto (\langle u_1, w\rangle, \ldots, \langle u_{k-1}, w\rangle)$ は、次元の比較から核に 0 でない元 $w$ をもちます。この $w$ は $c_1 = \cdots = c_{k-1} = 0$ なので $a(w,w) = \sum_{m\ge k}\lambda_mc_m^2 \ge \lambda_k\|w\|^2$、すなわち $\max_FR \ge \lambda_k$ です。

∎

これは、下に有界な自己共役作用素のミニマックス原理（定理 6.3［作用素の収束と摂動］）と、計数関数についての変分原理（補題 2.2［Dirichlet ラプラシアンの固有値と Weyl 則］）の、いまの設定での形です。行列の場合は Courant–Fischer の定理（定理 3.1［行列の固有値評価と不等式］）です。

定義 2.2（固有値問題の Galerkin 近似）↗

$V_h \subset V$ を $n$ 次元の部分空間とする。$u^h \in V_h$、$u^h \ne 0$ と実数 $\lambda^h$ が

$$a(u^h, v) = \lambda^h\langle u^h, v\rangle\qquad(\text{すべての}\ v \in V_h)$$

を満たすとき、$\lambda^h$ を Ritz 値、$u^h$ を Ritz ベクトルという。$V_h$ の基底 $\varphi_1, \ldots, \varphi_n$ について、剛性行列 $K = (a(\varphi_j, \varphi_i))$ と質量行列 $G = (\langle\varphi_j, \varphi_i\rangle)$ をとると、$u^h = \sum_jx_j\varphi_j$ について上の式は一般化固有値問題 $K\boldsymbol{x} = \lambda^hG\boldsymbol{x}$ である。$V_h \subset V$ であるとき近似は適合であるといい、そうでないとき（差分法、質量行列の集中化、不連続な試験関数など）非適合であるという。$a(w - P_hw, v) = 0$（すべての $v \in V_h$）で定まる $P_h : V \to V_h$ を Ritz 射影（$a$ についての直交射影）という。

Ritz 射影は 有限要素法 Iの 定義 2.2［有限要素法 I］ の $R_h$ と同じものです。$a$ についての直交射影なので、$\|w\|_a^2 = \|P_hw\|_a^2 + \|w - P_hw\|_a^2$、$\|w - P_hw\|_a = \min_{v\in V_h}\|w - v\|_a$ が成り立ちます（定理 2.5［有限要素法 I］ の対称な場合）。

補題 2.3（Ritz 値のミニマックス表示）✓

$V_h \subset V$ を $n$ 次元の部分空間とする。実数 $\lambda_1^h \le \cdots \le \lambda_n^h$ と $V_h$ の基底 $u_1^h, \ldots, u_n^h$ で

$$\langle u_i^h, u_j^h\rangle = \delta_{ij},\qquad a(u_i^h, u_j^h) = \lambda_i^h\delta_{ij}$$

を満たすものが存在し、Ritz 値はこの $\lambda_1^h, \ldots, \lambda_n^h$ ですべてである。さらに $k = 1, \ldots, n$ について

$$\lambda_k^h = \min_{\substack{F \subset V_h\\ \dim F = k}}\ \max_{w \in F,\ w \ne 0}R(w)$$

が成り立つ。

証明（補題 2.3）

$V_h$ の $\langle\cdot,\cdot\rangle$ についての正規直交基底 $\psi_1, \ldots, \psi_n$ をとり（定理 6.1［内積空間とグラム・シュミット直交化］）、$\hat K = (a(\psi_j, \psi_i))$ とします。$\hat K$ は実対称行列なので、直交行列 $O$ と対角行列 $\mathrm{diag}(\lambda_1^h, \ldots, \lambda_n^h)$（小さい順に並べる）で $O^{\mathsf T}\hat KO = \mathrm{diag}(\lambda_i^h)$ と書けます（系 4.3［スペクトル定理］）。$u_i^h = \sum_jO_{ji}\psi_j$ とおくと、$\langle u_i^h, u_j^h\rangle = (O^{\mathsf T}O)_{ij} = \delta_{ij}$、$a(u_i^h, u_j^h) = (O^{\mathsf T}\hat KO)_{ij} = \lambda_i^h\delta_{ij}$ です。$v = u_j^h$ で確かめると $a(u_i^h, v) = \lambda_i^h\langle u_i^h, v\rangle$ がすべての $v \in V_h$ で成り立つので、$\lambda_i^h$ は Ritz 値です。逆に $(\lambda^h, u^h)$ が Ritz 値と Ritz ベクトルなら、$u^h = \sum_ib_iu_i^h$ と書いて $v = u_j^h$ をとると $\lambda_j^hb_j = \lambda^hb_j$ で、$b_j \ne 0$ となる $j$ について $\lambda^h = \lambda_j^h$ です。

ミニマックス表示は、$w = \sum_ib_iu_i^h$ について $\|w\|^2 = \sum_ib_i^2$、$a(w,w) = \sum_i\lambda_i^hb_i^2$ となるので、補題 2.1 の証明が $u_m$ を $u_i^h$ に、$V$ を $V_h$ に替えてそのまま通ります。

∎

3. 差分法の固有値

定理 3.1（3 点差分の固有値）✓

$M \ge 2$、$h = 1/M$ とし、区間のモデル問題を 3 点差分で離散化した $(M-1)$ 次の行列 $A_h = \frac{1}{h^2}\mathrm{tridiag}(-1, 2, -1)$ を考える。

1. $A_h$ の固有値は $\lambda_k^{\mathrm{FD}} = \dfrac{4}{h^2}\sin^2\dfrac{k\pi h}{2}$（$k = 1, \ldots, M-1$）で、固有ベクトルは厳密な固有関数の格子値 $(\sin k\pi x_j)_{j=1}^{M-1}$ である。
2. $k = 1, \ldots, M-1$ について $$0 \lt \lambda_k - \lambda_k^{\mathrm{FD}} \le \frac{\lambda_k^2h^2}{12},\qquad \frac{4}{\pi^2} \lt \frac{\lambda_k^{\mathrm{FD}}}{\lambda_k} \lt 1.$$

証明（定理 3.1）

1. 補題 8.3［行列のノルムと線形方程式の条件数］ を $n = M - 1$、$\alpha = 2/h^2$、$\beta = -1/h^2$ で使うと、固有値は $\frac{2}{h^2}(1 - \cos\frac{k\pi}{M}) = \frac{4}{h^2}\sin^2\frac{k\pi h}{2}$、固有ベクトルは $(\sin(jk\pi/M))_j = (\sin k\pi x_j)_j$ です（命題 3.1［楕円型・放物型方程式の差分法］ の 1 と同じ内容です）。

2. $x = k\pi h/2$ とおくと $0 \lt x \lt \pi/2$ で、$\lambda_k = 4x^2/h^2$、$\lambda_k^{\mathrm{FD}} = 4\sin^2x/h^2$ です。$0 \lt \sin x \lt x$ から $\lambda_k^{\mathrm{FD}} \lt \lambda_k$。$x \gt 0$ で $\sin x \ge x - x^3/6$（Taylor の定理 定理 5.3［平均値の定理とテイラーの定理］ で、剰余項 $\frac{\cos\xi}{120}x^5$ は $0 \lt x \lt \pi/2$ で正）で、右辺は $x \lt \sqrt6$ で正なので、$\sin^2x \ge x^2 - x^4/3 + x^6/36 \ge x^2 - x^4/3$ です。よって

$$\lambda_k - \lambda_k^{\mathrm{FD}} = \frac{4}{h^2}(x^2 - \sin^2x) \le \frac{4}{h^2}\cdot\frac{x^4}{3} = \frac{k^4\pi^4h^2}{12} = \frac{\lambda_k^2h^2}{12}$$

です。比は $(\sin x/x)^2$ で、$\sin x/x$ は $(0, \pi/2]$ で減少し $x = \pi/2$ で $2/\pi$ なので、$4/\pi^2$ より大きくなります。

∎

差分の固有ベクトルが厳密な固有関数の格子値に一致するのは、定数係数で格子が等間隔という特別な事情によります。固有値の誤差の式から、相対誤差は $\lambda_kh^2/12 = (k\pi h)^2/12$ 以下です。相対誤差を 1% 以下にするには $k\pi h \le 0.35$、すなわち $k \le 0.11\,M$ 程度でなければなりません。$M - 1$ 個の固有値のうち、使えるのは下から 1 割ほどです。最大の固有値は $\lambda_{M-1}^{\mathrm{FD}} \approx 4/h^2$ で、対応する厳密な値 $\pi^2/h^2$ の約 4 割しかありません。

命題 3.2（5 点差分の固有値）✓

単位正方形の Dirichlet ラプラシアン（固有値 $\lambda_{jk} = \pi^2(j^2 + k^2)$、$j, k \ge 1$）の、$h = 1/M$ の格子での 5 点差分の固有値は $\lambda_{jk}^{\mathrm{FD}} = \lambda_j^{\mathrm{FD}} + \lambda_k^{\mathrm{FD}}$（$1 \le j, k \le M-1$）で、

$$0 \lt \lambda_{jk} - \lambda_{jk}^{\mathrm{FD}} \le \frac{\pi^4(j^4 + k^4)h^2}{12} \le \frac{\lambda_{jk}^2h^2}{12}$$

を満たす。両方を小さい順に並べたとき、$m$ 番目の差分の固有値は $m$ 番目の厳密な固有値より小さい（$m \le (M-1)^2$）。

証明（命題 3.2）

固有値が 1 次元の固有値の和であることは 命題 3.1［楕円型・放物型方程式の差分法］ の 2 です（固有ベクトルは $\sin(j\pi x_a)\sin(k\pi y_b)$）。厳密な固有値は 補題 3.2［Dirichlet ラプラシアンの固有値と Weyl 則］ によります。定理 3.1 の 2 を $j$ と $k$ について足すと不等式が出ます（$j^4 + k^4 \le (j^2 + k^2)^2$）。

並べ替えについて: 厳密な固有値を小さい順に $\mu_1 \le \mu_2 \le \cdots$ とし、$\mu_1, \ldots, \mu_m$ を与える添字 $(j,k)$ を $m$ 個選びます。どれも $\lambda_{jk} \le \mu_m$ を満たします。(a) 選んだ添字がすべて $1 \le j, k \le M-1$ の範囲にあるとき、それぞれについて $\lambda_{jk}^{\mathrm{FD}} \lt \lambda_{jk} \le \mu_m$ なので、$\mu_m$ より小さい差分の固有値が $m$ 個以上あり、$m$ 番目の差分の固有値は $\mu_m$ より小さくなります。(b) 範囲の外の添字（$j \ge M$ または $k \ge M$）が含まれるとき、$\mu_m \ge \pi^2(M^2 + 1)$ です。一方、差分の固有値はすべて $\frac{4}{h^2} + \frac{4}{h^2} = 8M^2$ より小さく、$8 \lt \pi^2$ なので $8M^2 \lt \mu_m$ です。

∎

「差分の固有値は下界になる」は、定数係数のラプラシアンと長方形の格子に特有の性質です。一般には成り立ちません。

例 3.3（差分の固有値が上にずれる例）

$-((1 + 10x)u')' = \lambda u$、$u(0) = u(1) = 0$ を、係数を格子点の中点で評価する標準的な 3 点差分

$$-\frac{p_{j+1/2}(u_{j+1} - u_j) - p_{j-1/2}(u_j - u_{j-1})}{h^2} = \lambda u_j,\qquad p_{j\pm1/2} = 1 + 10\,(x_j \pm h/2)$$

で離散化します。参照値は Chebyshev 選点法（スペクトル法の 命題 7.5［スペクトル法］ の微分行列）で求め、点の数を変えて一致を確かめます。

# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1
import numpy as np
from scipy.linalg import eigh_tridiagonal

def cheb(N):
    j = np.arange(N + 1)
    x = np.cos(np.pi * j / N)
    c = np.where((j == 0) | (j == N), 2.0, 1.0) * (-1.0) ** j
    D = (c[:, None] / c[None, :]) / (x[:, None] - x[None, :] + np.eye(N + 1))
    D -= np.diag(D.sum(axis=1))
    return D, x

p = lambda x: 1 + 10 * x
refs = []
for N in [40, 60, 80]:                               # 参照値: Chebyshev 選点法（N を変えて一致を確かめる）
    D, t = cheb(N)
    y = (t + 1) / 2
    refs.append(np.sort(np.linalg.eigvals((-4 * D @ np.diag(p(y)) @ D)[1:N, 1:N]).real)[:2])
    print(f"Chebyshev N = {N}: lambda_1 = {refs[-1][0]:.10f}, lambda_2 = {refs[-1][1]:.10f}")
print("    M    lambda_1^h - lambda_1    lambda_2^h - lambda_2")
for M in [8, 16, 32, 64]:
    h = 1.0 / M
    pm = p(h * (np.arange(M) + 0.5))                 # 係数は格子点の中点で評価する
    fd = eigh_tridiagonal((pm[:-1] + pm[1:]) / h**2, -pm[1:-1] / h**2, eigvals_only=True)[:2]
    print(f"{M:5d}    {fd[0] - refs[-1][0]:+.4e}              {fd[1] - refs[-1][1]:+.4e}")

Chebyshev N = 40: lambda_1 = 44.4515837733, lambda_2 = 182.1541741516
Chebyshev N = 60: lambda_1 = 44.4515837733, lambda_2 = 182.1541741516
Chebyshev N = 80: lambda_1 = 44.4515837733, lambda_2 = 182.1541741516
    M    lambda_1^h - lambda_1    lambda_2^h - lambda_2
    8    +3.7299e-01              -6.8046e+00
   16    +1.0638e-01              -1.9742e+00
   32    +2.7831e-02              -5.1122e-01
   64    +7.0477e-03              -1.2891e-01

第 1 固有値は上から、第 2 固有値は下から近づいています。どちらも誤差は $h$ を半分にするごとに約 $1/4$（$O(h^2)$）ですが、向きは固有値ごとに違います。差分法は非適合な近似で、Rayleigh–Ritz の上界のような一方向の保証がありません。

例 3.4（離散固有値の個数と Weyl の法則）

単位正方形では、$\lambda$ 以下の固有値の個数 $N(\lambda)$ は $\lambda \to \infty$ で $\lambda/(4\pi)$ のように増えます（定理 5.2［Dirichlet ラプラシアンの固有値と Weyl 則］ で $n = 2$、$|\Omega| = 1$、$\omega_2 = \pi$）。$h = 1/64$ の 5 点差分は $63^2 = 3969$ 個の固有値をもち、すべて $8/h^2 = 32768$ より小さい値です。$\lambda$ 以下の個数を数えます。

# Python 3.9.6 / NumPy 2.0.2
import numpy as np

M = 64
h = 1.0 / M
k = np.arange(1, M)
lam1 = 4 * np.sin(k * np.pi * h / 2) ** 2 / h**2
fd = np.add.outer(lam1, lam1).ravel()                       # 5 点差分の固有値（(M-1)^2 個）
kk = np.arange(1, 400)
exact = np.add.outer((kk * np.pi) ** 2, (kk * np.pi) ** 2).ravel()
print(f"number of discrete eigenvalues = {fd.size}, largest = {fd.max():.1f}, 8/h^2 = {8 / h**2:.1f}")
print("  lambda    lambda h^2    exact N(lambda)    5-point N_h(lambda)    Weyl lambda/(4 pi)")
for lam in [500, 2000, 8000, 16000, 32000]:
    print(f"  {lam:6d}    {lam * h**2:8.3f}     {np.sum(exact <= lam):8d}           {np.sum(fd <= lam):8d}"
          f"             {lam / (4 * np.pi):10.1f}")

number of discrete eigenvalues = 3969, largest = 32748.3, 8/h^2 = 32768.0
  lambda    lambda h^2    exact N(lambda)    5-point N_h(lambda)    Weyl lambda/(4 pi)
     500       0.122           33                 33                   39.8
    2000       0.488          146                150                  159.2
    8000       1.953          612                705                  636.6
   16000       3.906         1234               1869                 1273.2
   32000       7.812         2490               3917                 2546.5

差分の固有値は真の値より小さいので（命題 3.2）、個数は多めに数えられます。$\lambda h^2 \approx 0.5$ までは差が数個ですが、$\lambda h^2 \approx 2$ で 15% 多く、それより上では離散固有値が $8/h^2$ の下に押し込められて、個数が真の値の 1.5 倍になります。作用素の固有値は無限個あり、行列の固有値は有限個です。上位の離散固有値は、作用素のスペクトルではなく格子の性質を映しています。

4. Ritz 値は上から近づく

定理 4.1（適合な Galerkin 近似の上界と単調性）✓

$V_h \subset V$ を $n$ 次元の部分空間、$\lambda_1^h \le \cdots \le \lambda_n^h$ をその Ritz 値とする。

1. $k = 1, \ldots, n$ について $\lambda_k \le \lambda_k^h$。
2. $V_h \subset V_{h'} \subset V$ で、$V_{h'}$ の Ritz 値を $\lambda_k^{h'}$ とすると、$k = 1, \ldots, n$ について $\lambda_k^{h'} \le \lambda_k^h$。

証明（定理 4.1）

補題 2.3 と 補題 2.1 から、$\lambda_k^h$、$\lambda_k^{h'}$、$\lambda_k$ は、同じ量 $\max_{w\in F}R(w)$ の、$k$ 次元部分空間 $F$ をそれぞれ $V_h$、$V_{h'}$、$V$ の中で動かしたときの最小値です。$V_h \subset V_{h'} \subset V$ なので、動かす範囲が広いほど最小値は小さく、$\lambda_k \le \lambda_k^{h'} \le \lambda_k^h$ です。

∎

1 は 定理 3.1［無限次元の固有値の包含］、系 6.4［作用素の収束と摂動］ と同じ内容で、2 は領域の単調性（系 3.1［Dirichlet ラプラシアンの固有値と Weyl 則］）やミニマックス値の単調性（系 6.5［作用素の収束と摂動］）と同じ仕組みです。有限要素法では、分割を細分して得た空間はもとの空間を含むので、細分するたびに各 Ritz 値は減ります（減らないこともあるが増えない）。行列の Lanczos 法で、Krylov 部分空間を広げると Ritz 値が単調に端の固有値へ近づくこと（定理 4.3［Krylov 部分空間法 I］）も同じ原理です。

次に、上からどれだけ近いかを評価します。主張は、固有値の誤差が固有関数の近似誤差の 2 乗になる、というものです。たとえば固有関数の誤差が 10 分の 1 になれば、固有値の誤差は 100 分の 1 になります。証明では、固有関数 $u$ の誤差を直接見るかわりに $A^{-1}u$ を挟み、Galerkin 直交性（誤差が $V_h$ と $a$ について直交すること）を使って、誤差の積の形にします。

定理 4.2（Ritz 値の誤差は固有関数のエネルギー誤差の 2 乗）✓

$V_h \subset V$ を有限次元の部分空間、$P_h$ を Ritz 射影、$k \ge 1$、$E_k = \mathrm{span}\{u_1, \ldots, u_k\}$ とする。

1. $\varepsilon_k = \max\{\|u - P_hu\|_a : u \in E_k,\ \|u\| = 1\}$ とおく。$2\varepsilon_k^2 \lt \lambda_1$ ならば $\dim V_h \ge k$ で $$\lambda_k \le \lambda_k^h \le \frac{\lambda_k}{1 - 2\varepsilon_k^2/\lambda_1}.$$
2. 定数 $\eta \gt 0$ があって、すべての $w \in D(A)$ について $\|w - P_hw\|_a \le \eta\|Aw\|$ が成り立つとする。$2\eta^2\lambda_k \lt 1$ ならば $\dim V_h \ge k$ で $$\lambda_k \le \lambda_k^h \le \frac{\lambda_k}{1 - 2\eta^2\lambda_k},\qquad\text{とくに}\ 4\eta^2\lambda_k \le 1\ \text{なら}\quad 0 \le \lambda_k^h - \lambda_k \le 4\eta^2\lambda_k^2.$$

証明（定理 4.2）

$u = \sum_{m\le k}c_mu_m \in E_k$、$\|u\| = 1$ とし、$z = A^{-1}u = \sum_{m\le k}(c_m/\lambda_m)u_m \in E_k$ とおきます。

内積 $\langle u, u - P_hu\rangle$ の評価。 $Az = u$ なので $\langle u, u - P_hu\rangle = a(z, u - P_hu)$ です。Ritz 射影の定義から $a(P_hz, u - P_hu) = 0$ なので、Cauchy–Schwarz の不等式（$a$ は $V$ の内積）から

$$\langle u, u - P_hu\rangle = a(z - P_hz, u - P_hu) \le \|z - P_hz\|_a\,\|u - P_hu\|_a$$

です。1 の場合: $z \in E_k$、$\|z\|^2 = \sum c_m^2/\lambda_m^2 \le 1/\lambda_1^2$ なので $\|z - P_hz\|_a \le \varepsilon_k\|z\| \le \varepsilon_k/\lambda_1$、$\|u - P_hu\|_a \le \varepsilon_k$ で、右辺は $\sigma := \varepsilon_k^2/\lambda_1$ 以下です。2 の場合: $\|z - P_hz\|_a \le \eta\|Az\| = \eta$、$\|u - P_hu\|_a \le \eta\|Au\| \le \eta\lambda_k$（$\|Au\|^2 = \sum_{m\le k}\lambda_m^2c_m^2 \le \lambda_k^2$）で、右辺は $\sigma := \eta^2\lambda_k$ 以下です。

$P_hu$ の大きさ。 どちらの場合も

$$\|P_hu\|^2 = \|u - (u - P_hu)\|^2 = 1 - 2\langle u, u - P_hu\rangle + \|u - P_hu\|^2 \ge 1 - 2\sigma \gt 0$$

です。したがって $P_h$ は $E_k$ の上で単射で、$F = P_hE_k \subset V_h$ は $k$ 次元です（とくに $\dim V_h \ge k$）。

Rayleigh 商。 $P_h$ は $a$ についての直交射影なので $a(P_hu, P_hu) \le a(u, u) \le \lambda_k$ です。よって $F$ の 0 でない元 $P_hu$（$u \in E_k$ を $\|u\| = 1$ に正規化）について $R(P_hu) \le \lambda_k/(1 - 2\sigma)$ で、補題 2.3 から $\lambda_k^h \le \max_FR \le \lambda_k/(1 - 2\sigma)$ です。下界 $\lambda_k \le \lambda_k^h$ は 定理 4.1 です。$4\eta^2\lambda_k \le 1$ なら $1 - 2\sigma \ge 1/2$ で、$\lambda_k^h - \lambda_k \le \lambda_k\cdot2\sigma/(1 - 2\sigma) \le 4\sigma\lambda_k = 4\eta^2\lambda_k^2$ です。

∎

証明の要点は、$\langle u, u - P_hu\rangle$ が誤差の 2 つの因子の積になることです。$u - P_hu$ を素直に評価すると $\|u - P_hu\| = O(\varepsilon)$ ですが、Galerkin 直交性のおかげで $z - P_hz$ をもう 1 つ挟め、$O(\varepsilon^2)$ になります。境界値問題の $L^2$ 誤差の評価（Aubin–Nitsche の双対論法。定理 5.2［有限要素法 I］）と同じ仕組みです。Rayleigh 商が固有ベクトルで停留する（1 次の変分が 0）ことの現れでもあります。

2 の仮定は、有限要素法では補間誤差の評価と $H^2$ 正則性（ラプラシアンの値で 2 階微分の大きさが押さえられること）から出ます。

• 1 次元の区分 1 次要素（幅の最大が $h$ の任意の分割）: Ritz 射影はエネルギーノルムでの最良近似なので、節点補間 $I_hw \in V_h$ と比べて、$w \in D(A) = H^2\cap H^1_0$ について $\|w - P_hw\|_a \le \|(w - I_hw)'\| \le \frac h\pi\|w''\|$（定理 4.1［有限要素法 I］ の 1）で、$\eta = h/\pi$ ととれます（等間隔の分割では Ritz 射影は節点補間そのものです。補題 6.2［無限次元の固有値の包含］）。したがって $\lambda_kh^2 \le \pi^2/4$（$kh \le 1/2$）のとき $$0 \le \lambda_k^h - \lambda_k \le \frac{4}{\pi^2}\lambda_k^2h^2.$$
• 2 次元の三角形要素: 単位正方形（$C_{\mathrm{reg}} = 1$。より一般に、$H^2$ 正則性 $\|D^2w\| \le C_{\mathrm{reg}}\|\Delta w\|$ が成り立つ多角形領域）の形状正則な分割では、Céa の補題と補間誤差の評価から $\|w - P_hw\|_a \le C_Ih\,|w|_{H^2} \le C_IC_{\mathrm{reg}}h\,\|Aw\|$ で（有限要素法 IIの 定理 5.1［有限要素法 II］ と 補題 5.2［有限要素法 II］）、$\eta = C_IC_{\mathrm{reg}}h$ です。固有値の誤差はやはり $O(\lambda_k^2h^2)$ です。

どちらも、誤差の定数が $\lambda_k^2$ に比例します。差分法の場合（定理 3.1）と同じで、高い固有値ほど細かい格子が要ります。固有関数が $H^2$ に属さない領域（凹角のある多角形）では 2 の仮定が成り立たず、1 の $\varepsilon_k$ が $O(h)$ より遅くしか減らないので、次数が落ちます（例 4.5）。

1 次元の等間隔の分割では、Ritz 値が閉じた形で求まり、誤差の定数まで分かります。

定理 4.3（1 次元の区分 1 次要素の固有値）✓

区間のモデル問題を、等間隔の分割（$h = 1/M$）の区分 1 次の有限要素空間 $V_h \subset H^1_0(0,1)$ で離散化する。$\theta_k = k\pi h$（$k = 1, \ldots, M-1$）とおく。

1. Ritz 値は $$\lambda_k^h = \frac{6}{h^2}\cdot\frac{1 - \cos\theta_k}{2 + \cos\theta_k}\qquad(k = 1, \ldots, M-1)$$ で、Ritz ベクトルは $\sin k\pi x$ の節点補間である。
2. 質量行列 $G$ を対角行列 $hI$ に置き換える（集中化）と、一般化固有値問題は 3 点差分の固有値問題に一致し、固有値は 定理 3.1 の $\lambda_k^{\mathrm{FD}}$ になる。
3. $k$ を固定して $h \to 0$ とすると $$\lambda_k^h = \lambda_k + \frac{\lambda_k^2h^2}{12} + \frac{\lambda_k^3h^4}{360} + O(h^6),\qquad \lambda_k^{\mathrm{FD}} = \lambda_k - \frac{\lambda_k^2h^2}{12} + \frac{\lambda_k^3h^4}{360} + O(h^6)$$ で、平均は $\frac12(\lambda_k^h + \lambda_k^{\mathrm{FD}}) = \lambda_k + \dfrac{\lambda_k^3h^4}{360} + O(h^6)$ である。

証明（定理 4.3）

1. 節点基底についての剛性行列と質量行列は $K = \frac1h\mathrm{tridiag}(-1,2,-1)$、$G = \frac h6\mathrm{tridiag}(1,4,1)$ です（命題 6.1［有限要素法 I］ の 1）。補題 8.3［行列のノルムと線形方程式の条件数］ から、$\boldsymbol{v}_k = (\sin j\theta_k)_{j=1}^{M-1}$ は $K$ と $G$ の共通の固有ベクトルで

$$K\boldsymbol{v}_k = \frac{2 - 2\cos\theta_k}{h}\boldsymbol{v}_k,\qquad G\boldsymbol{v}_k = \frac{h(4 + 2\cos\theta_k)}{6}\boldsymbol{v}_k$$

です。よって $K\boldsymbol{v}_k = \lambda_k^hG\boldsymbol{v}_k$、$\lambda_k^h = \frac{6(2 - 2\cos\theta_k)}{h^2(4 + 2\cos\theta_k)}$ です。$\boldsymbol{v}_1, \ldots, \boldsymbol{v}_{M-1}$ は 1 次独立で、$V_h$ の次元は $M - 1$ なので、Ritz 値はこれですべてです（補題 2.3）。$\boldsymbol{v}_k$ を係数とする $V_h$ の元は、節点 $x_j$ で $\sin k\pi x_j$ をとる区分 1 次関数です。

2. $K\boldsymbol{x} = \lambda hI\boldsymbol{x}$ は $\frac{1}{h^2}\mathrm{tridiag}(-1,2,-1)\boldsymbol{x} = \lambda\boldsymbol{x}$ です。

3. $\theta = \theta_k$ とすると $\theta^2 = \lambda_kh^2$ で、$\lambda_k^{\mathrm{FD}}/\lambda_k = \Psi(\theta) := 2(1 - \cos\theta)/\theta^2$、$\lambda_k^h/\lambda_k = \Phi(\theta) := 3\Psi(\theta)/(2 + \cos\theta)$ です。Taylor の定理（定理 5.3［平均値の定理とテイラーの定理］）から $1 - \cos\theta = \frac{\theta^2}{2} - \frac{\theta^4}{24} + \frac{\theta^6}{720} + O(\theta^8)$ なので

$$\Psi(\theta) = 1 - \frac{\theta^2}{12} + \frac{\theta^4}{360} + O(\theta^6),\qquad \frac{3}{2 + \cos\theta} = \Bigl(1 - \frac{\theta^2}{6} + \frac{\theta^4}{72} + O(\theta^6)\Bigr)^{-1} = 1 + \frac{\theta^2}{6} + \frac{\theta^4}{72} + O(\theta^6)$$

です（最後は $(1 - c)^{-1} = 1 + c + c^2 + O(c^3)$ に $c = \theta^2/6 - \theta^4/72$ を入れ、$\theta^4$ の係数 $-\frac1{72} + \frac1{36} = \frac1{72}$ を得る）。積をとると

$$\Phi(\theta) = 1 + \Bigl(\frac16 - \frac1{12}\Bigr)\theta^2 + \Bigl(\frac1{72} - \frac1{72} + \frac1{360}\Bigr)\theta^4 + O(\theta^6) = 1 + \frac{\theta^2}{12} + \frac{\theta^4}{360} + O(\theta^6)$$

です。$\lambda_k$ を掛け、$\theta^2 = \lambda_kh^2$ を入れると主張の式になります。

∎

適合な近似（一致質量）は上から、集中質量（＝差分）は下から、同じ大きさ $\lambda_k^2h^2/12$ の誤差で近づきます。先頭の誤差項が打ち消し合うので、2 つの平均は 4 次の精度になります。定理 4.2 から出る上界の定数 $4/\pi^2 \approx 0.41$ は、正確な定数 $1/12 \approx 0.083$ の約 5 倍で、次数 $\lambda_k^2h^2$ は一致しています。

例 4.4（1 次元の 3 つの離散化の比較）

# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1
import numpy as np
from scipy.linalg import eigh, eigh_tridiagonal

def eigs(M):
    """-u'' = lambda u（u(0) = u(1) = 0）、h = 1/M。差分・一致質量の有限要素・集中質量の固有値。"""
    h, n = 1.0 / M, M - 1
    K = (2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)) / h
    Mass = h * (4 * np.eye(n) + np.eye(n, k=1) + np.eye(n, k=-1)) / 6
    fd = eigh_tridiagonal(2 * np.ones(n) / h**2, -np.ones(n - 1) / h**2, eigvals_only=True)
    fem = eigh(K, Mass, eigvals_only=True)
    lumped = eigh(K, h * np.eye(n), eigvals_only=True)
    return fd, fem, lumped

k = np.arange(1, 6)
exact = (k * np.pi) ** 2
for M in [16, 64]:
    h = 1.0 / M
    fd, fem, lumped = eigs(M)
    th = k * np.pi * h
    closed_fd = 4 * np.sin(th / 2) ** 2 / h**2
    closed_fem = 6 * (1 - np.cos(th)) / (h**2 * (2 + np.cos(th)))
    print(f"M = {M}: max |closed form - computed|: FD {np.max(np.abs(fd[:5] - closed_fd)):.1e}, "
          f"FEM {np.max(np.abs(fem[:5] - closed_fem)):.1e}, lumped - FD {np.max(np.abs(lumped[:5] - fd[:5])):.1e}")
    print("   k    FD - exact     FEM - exact    average - exact   (FEM - exact) / (lambda^2 h^2 / 12)")
    for i in range(5):
        print(f"   {k[i]}    {fd[i] - exact[i]:+.4e}    {fem[i] - exact[i]:+.4e}    "
              f"{0.5 * (fd[i] + fem[i]) - exact[i]:+.4e}       {(fem[i] - exact[i]) / (exact[i]**2 * h**2 / 12):.4f}")
fd, fem, _ = eigs(64)
print("M = 64, k = 63 (the largest): FD / exact =", f"{fd[-1] / (63 * np.pi)**2:.4f}",
      ", FEM / exact =", f"{fem[-1] / (63 * np.pi)**2:.4f}")

M = 16: max |closed form - computed|: FD 8.8e-13, FEM 2.6e-13, lumped - FD 9.8e-13
   k    FD - exact     FEM - exact    average - exact   (FEM - exact) / (lambda^2 h^2 / 12)
   1    -3.1668e-02    +3.1749e-02    +4.0655e-05       1.0013
   2    -5.0474e-01    +5.0991e-01    +2.5834e-03       1.0051
   3    -2.5389e+00    +2.5970e+00    +2.9057e-02       1.0111
   4    -7.9523e+00    +8.2726e+00    +1.6013e-01       1.0191
   5    -1.9192e+01    +2.0380e+01    +5.9395e-01       1.0284
M = 64: max |closed form - computed|: FD 4.2e-12, FEM 6.2e-12, lumped - FD 4.7e-12
   k    FD - exact     FEM - exact    average - exact   (FEM - exact) / (lambda^2 h^2 / 12)
   1    -1.9816e-03    +1.9820e-03    +1.5916e-07       1.0001
   2    -3.1699e-02    +3.1719e-02    +1.0181e-05       1.0003
   3    -1.6041e-01    +1.6064e-01    +1.1589e-04       1.0007
   4    -5.0669e-01    +5.0799e-01    +6.5048e-04       1.0013
   5    -1.2361e+00    +1.2411e+00    +2.4781e-03       1.0020
M = 64, k = 63 (the largest): FD / exact = 0.4180 , FEM / exact = 1.2525

符号は定理のとおりで、$h$ を $1/4$ にすると差分と有限要素の誤差は約 $1/16$、平均の誤差は約 $1/256$ になっています（$k = 1$ で $4.07\times10^{-5} \to 1.59\times10^{-7}$。予言値 $\lambda_1^3h^4/360$ は $4.07\times10^{-5}$ と $1.59\times10^{-7}$）。誤差は $k$ とともに $\lambda_k^2 \propto k^4$ で増え、$M = 16$ の第 5 固有値（$\lambda_5 \approx 247$）は 8% ずれています。最大の離散固有値は、差分では厳密な値の $0.42$ 倍、有限要素では $1.25$ 倍です。

例 4.5（単位正方形と L 字形領域）

単位正方形の固有値は $\pi^2(j^2 + k^2)$ で、小さい方から $2\pi^2$、$5\pi^2$（2 重）、$8\pi^2$、$10\pi^2$（2 重）です。正方形の格子を対角線で直角二等辺三角形に分けた区分 1 次要素で、下位 6 個の Ritz 値を求めます（scipy.sparse.linalg.eigsh のシフト反転、$\sigma = 0$）。あわせて、1 辺 1 の正方形 3 個からなる L 字形領域の第 1 固有値を求めます。参照値 $9.6397238440219$ は Betcke–Trefethen の計算（参考文献 [5]）によります。

# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import eigsh

def p1_matrices(cells, h):
    """cells[i, j] が真の正方形セル（1 辺 h）を対角線で 2 つの直角二等辺三角形に分け、
    区分 1 次要素の剛性行列 K と質量行列 M を内部節点について組み立てる。"""
    n = cells.shape[0]
    i, j = np.nonzero(cells)
    idx = lambda a, b: a * (n + 1) + b
    tris = np.concatenate([np.stack([idx(i, j), idx(i + 1, j), idx(i + 1, j + 1)], axis=1),
                           np.stack([idx(i, j), idx(i, j + 1), idx(i + 1, j + 1)], axis=1)])
    KE = np.array([[1, -1, 0], [-1, 2, -1], [0, -1, 1]]) / 2.0     # 直角は 2 番目の頂点
    ME = h**2 * (np.ones((3, 3)) + np.eye(3)) / 24.0
    rows, cols = np.repeat(tris, 3, axis=1).ravel(), np.tile(tris, (1, 3)).ravel()
    size = (n + 1) ** 2
    K = coo_matrix((np.tile(KE.ravel(), len(tris)), (rows, cols)), shape=(size, size)).tocsr()
    M = coo_matrix((np.tile(ME.ravel(), len(tris)), (rows, cols)), shape=(size, size)).tocsr()
    inner = cells[:-1, :-1] & cells[:-1, 1:] & cells[1:, :-1] & cells[1:, 1:]   # まわりの 4 セルが領域内
    a, b = np.nonzero(inner)
    keep = idx(a + 1, b + 1)
    return K[keep][:, keep], M[keep][:, keep]

exact = np.pi**2 * np.array([2, 5, 5, 8, 10, 10])
print("unit square: (lambda_h - lambda) / (lambda^2 h^2) for the six smallest eigenvalues")
for n in [16, 32, 64]:
    h = 1.0 / n
    K, M = p1_matrices(np.ones((n, n), dtype=bool), h)
    lam = np.sort(eigsh(K, k=6, M=M, sigma=0.0, which="LM", return_eigenvectors=False))
    print(f"  h = 1/{n:<3d}: " + "  ".join(f"{(l - e) / (e**2 * h**2):.4f}" for l, e in zip(lam, exact)))
print("  h = 1/64 : lambda_h =", "  ".join(f"{l:.5f}" for l in lam))

ref = 9.6397238440219          # L 字形領域（1 辺 1 の正方形 3 個）の第 1 固有値の文献値
print("L-shaped domain: the smallest eigenvalue")
prev = None
for n in [16, 32, 64, 128]:
    cells = np.ones((n, n), dtype=bool)
    cells[n // 2:, n // 2:] = False
    K, M = p1_matrices(cells, 2.0 / n)
    err = eigsh(K, k=1, M=M, sigma=0.0, which="LM", return_eigenvectors=False)[0] - ref
    print(f"  h = 1/{n // 2:<3d}: lambda_h = {ref + err:.6f}, error = {err:.3e}"
          + ("" if prev is None else f", order = {np.log2(prev / err):.2f}"))
    prev = err

unit square: (lambda_h - lambda) / (lambda^2 h^2) for the six smallest eigenvalues
  h = 1/16 : 0.1252  0.0860  0.1351  0.1238  0.0989  0.1012
  h = 1/32 : 0.1251  0.0860  0.1343  0.1247  0.0985  0.0990
  h = 1/64 : 0.1250  0.0860  0.1341  0.1249  0.0984  0.0985
  h = 1/64 : lambda_h = 19.75110  49.39914  49.42774  79.14698  98.92999  98.93031
L-shaped domain: the smallest eigenvalue
  h = 1/8  : lambda_h = 9.965977, error = 3.263e-01
  h = 1/16 : lambda_h = 9.740817, error = 1.011e-01, order = 1.69
  h = 1/32 : lambda_h = 9.672951, error = 3.323e-02, order = 1.61
  h = 1/64 : lambda_h = 9.651203, error = 1.148e-02, order = 1.53

単位正方形では、誤差を $\lambda^2h^2$ で割った値が $h$ によらず $0.09$〜$0.13$ で、定理 4.2 の $O(\lambda_k^2h^2)$ のとおりです。誤差はすべて正です（定理 4.1）。厳密には 2 重の固有値 $5\pi^2 \approx 49.348$ が、離散では $49.399$ と $49.428$ に分かれています。対角線の向きが分割の対称性を壊し、2 つの固有関数（の 1 次結合）の近似の精度が異なるためです。

L 字形領域では、次数が 2 より小さく、$h$ とともに下がっています。凹角（内角 $3\pi/2$）の頂点の近くで固有関数が $r^{2/3}\sin(2\vartheta/3)$ のように振る舞い、$H^2$ に属さないためです。一様な分割ではエネルギー誤差が $O(h^{2/3})$、固有値の誤差がその 2 乗の $O(h^{4/3})$ になることが知られていて（Strang–Fix の教科書、参考文献 [1]）、観測された次数はこの値に向かって下がっています。頂点に向けて分割を細かくすれば次数 2 が回復します。

5. 固有ベクトルの誤差

固有値の誤差は固有値そのものだけで決まりましたが、固有ベクトルの誤差は、隣の固有値とのすき間に依存します。

次の命題の証明の方針は、Ritz 射影した $P_hu_k$ を「近似の固有ベクトル」とみて、Ritz ベクトルの基底で展開することです。$k$ 番目以外の成分は、鍵になる等式により（誤差）÷（固有値の差）で押さえられます。この「固有値の差」が小さいと割り算で誤差が大きくなる、というのが、すき間 $g_k$ が定数に入る理由です。

補題 5.1（Ritz 射影の誤差の双対評価）✓

定理 4.2 の 2 の仮定（すべての $w \in D(A)$ で $\|w - P_hw\|_a \le \eta\|Aw\|$）のもとで、すべての $w \in V$ について

$$\|w - P_hw\| \le \eta\,\|w - P_hw\|_a$$

が成り立つ。とくに $w \in D(A)$ なら $\|w - P_hw\| \le \eta^2\|Aw\|$ である。

証明（補題 5.1）

$e = w - P_hw$、$z = A^{-1}e \in D(A)$ とおくと、$\|e\|^2 = \langle Az, e\rangle = a(z, e)$ です。$a(P_hz, e) = a(e, P_hz) = 0$（Ritz 射影の定義で $v = P_hz \in V_h$）なので

$$\|e\|^2 = a(z - P_hz, e) \le \|z - P_hz\|_a\|e\|_a \le \eta\|Az\|\,\|e\|_a = \eta\|e\|\,\|e\|_a$$

で、$\|e\|$ で割ると主張が出ます（$e = 0$ のときは自明）。

∎

命題 5.2（単純固有値の Ritz ベクトルの誤差）✓

定理 4.2 の 2 の仮定のもとで、$n = \dim V_h \ge k$ とし、$\lambda_k$ を単純固有値（$\lambda_{k-1} \lt \lambda_k \lt \lambda_{k+1}$。$k = 1$ では右の不等式だけ）、$u_k$ を $\|u_k\| = 1$ の固有ベクトルとする。補題 2.3 の Ritz ベクトル $u_j^h$ について、$j \ne k$ で $\lambda_j^h \ne \lambda_k$ とし、

$$\rho = \max_{j \ne k}\frac{\lambda_k}{|\lambda_j^h - \lambda_k|}$$

とおく。$u_k^h$ の符号を $\langle P_hu_k, u_k^h\rangle \ge 0$ となるように選ぶと

$$\|u_k - u_k^h\| \le 2(1 + \rho)\,\|u_k - P_hu_k\| \le 2(1 + \rho)\,\eta^2\lambda_k,\qquad \|u_k - u_k^h\|_a^2 = \lambda_k\|u_k - u_k^h\|^2 + (\lambda_k^h - \lambda_k)$$

が成り立つ。さらに、すき間 $g_k = \min(\lambda_k - \lambda_{k-1}, \lambda_{k+1} - \lambda_k)$（$k = 1$ では $\lambda_2 - \lambda_1$）について $8\eta^2\lambda_k^2 \le g_k$ かつ $4\eta^2\lambda_k \le 1$ ならば、$\rho \le 2\lambda_k/g_k$ で

$$\|u_k - u_k^h\|_a \le 2\eta\lambda_k\Bigl(1 + \Bigl(1 + \frac{2\lambda_k}{g_k}\Bigr)\eta\lambda_k^{1/2}\Bigr)$$

である。

証明（命題 5.2）

鍵になる等式。 $j = 1, \ldots, n$ について、Ritz ベクトルの式（$v = P_hu_k \in V_h$）と $a$ の対称性、Ritz 射影の定義、$a(u_k, v) = \lambda_k\langle u_k, v\rangle$ から

$$\lambda_j^h\langle P_hu_k, u_j^h\rangle = a(P_hu_k, u_j^h) = a(u_k, u_j^h) = \lambda_k\langle u_k, u_j^h\rangle$$

です。両辺から $\lambda_k\langle P_hu_k, u_j^h\rangle$ を引くと $(\lambda_j^h - \lambda_k)\langle P_hu_k, u_j^h\rangle = \lambda_k\langle u_k - P_hu_k, u_j^h\rangle$ です。

$L^2$ ノルム。 $\beta_j = \langle P_hu_k, u_j^h\rangle$ とおくと $P_hu_k = \sum_j\beta_ju_j^h$（$\{u_j^h\}$ は $V_h$ の正規直交基底）で、上の等式と $\rho$ の定義、Bessel の不等式から

$$\|P_hu_k - \beta_ku_k^h\|^2 = \sum_{j\ne k}\beta_j^2 \le \rho^2\sum_{j\ne k}\langle u_k - P_hu_k, u_j^h\rangle^2 \le \rho^2\|u_k - P_hu_k\|^2$$

です。よって $\|u_k - \beta_ku_k^h\| \le (1 + \rho)\|u_k - P_hu_k\|$ です。$\beta_k \ge 0$ なので $|1 - \beta_k| = \bigl|\|u_k\| - \|\beta_ku_k^h\|\bigr| \le \|u_k - \beta_ku_k^h\|$ で、

$$\|u_k - u_k^h\| \le \|u_k - \beta_ku_k^h\| + |\beta_k - 1|\,\|u_k^h\| \le 2(1 + \rho)\|u_k - P_hu_k\|$$

です。補題 5.1 から $\|u_k - P_hu_k\| \le \eta^2\|Au_k\| = \eta^2\lambda_k$ です。

エネルギーノルム。 $a(u_k, u_k) = \lambda_k$、$a(u_k^h, u_k^h) = \lambda_k^h$、$a(u_k, u_k^h) = \lambda_k\langle u_k, u_k^h\rangle$、$\|u_k - u_k^h\|^2 = 2 - 2\langle u_k, u_k^h\rangle$ から

$$\|u_k - u_k^h\|_a^2 = \lambda_k - 2\lambda_k\langle u_k, u_k^h\rangle + \lambda_k^h = \lambda_k\|u_k - u_k^h\|^2 + (\lambda_k^h - \lambda_k)$$

です。

すき間による評価。 $j \gt k$ なら 定理 4.1 から $\lambda_j^h \ge \lambda_j \ge \lambda_{k+1}$ で、$\lambda_j^h - \lambda_k \ge g_k$。$j \lt k$ なら、$4\eta^2\lambda_{k-1} \le 4\eta^2\lambda_k \le 1$ なので 定理 4.2 の 2 から $\lambda_j^h \le \lambda_{k-1}^h \le \lambda_{k-1} + 4\eta^2\lambda_{k-1}^2 \le \lambda_{k-1} + g_k/2$ で、$\lambda_k - \lambda_j^h \ge g_k/2$。よって $\rho \le 2\lambda_k/g_k$ です。エネルギーノルムの等式に、$L^2$ の評価と $\lambda_k^h - \lambda_k \le 4\eta^2\lambda_k^2$ を入れ、$\sqrt{x^2 + y^2} \le x + y$ を使うと

$$\|u_k - u_k^h\|_a \le \lambda_k^{1/2}\cdot2\Bigl(1 + \frac{2\lambda_k}{g_k}\Bigr)\eta^2\lambda_k + 2\eta\lambda_k$$

で、主張の式です。

∎

区分 1 次の有限要素法（$\eta = O(h)$）では、Ritz ベクトルの誤差はエネルギーノルムで $O(h)$、$L^2$ ノルムで $O(h^2)$ です。境界値問題の誤差と同じ次数で、固有値の誤差 $O(h^2)$ はエネルギーノルムの誤差の 2 乗になっています（エネルギーノルムの等式がこの関係を厳密に表しています）。定数には $\lambda_k/g_k$ が入り、近接した固有値の固有ベクトルは精度が悪くなります。これは対称行列の固有ベクトルの摂動についての $\sin\theta$ 評価（定理 5.2［対称固有値問題と特異値分解の計算］。残差をすき間で割る）と同じ構造で、鍵になる等式は「$P_hu_k$ を近似固有ベクトルとみたときの残差」を表しています。例 4.5 で 2 重固有値が分かれたように、重複した固有値では個々の固有ベクトルは定まらず、固有空間どうしの距離を考える必要があります。

注意 5.3（重複固有値と固有値の集まり（文献））

重複した固有値、近接した固有値の集まり、非対称な問題を含む一般論は Babuška–Osborn の理論で、固有空間のあいだのすき間（gap）と、固有値の平均の誤差が評価されます（参考文献 [2]）。証明は作用素の Riesz 射影と解作用素の収束を使うので、この章では単純固有値の場合に限りました。

6. 上界と下界で挟む

Ritz 値は上界です。同じ仮定から下界も出るので、真の固有値を挟めます。

命題 6.1（Ritz 値による上界と Liu の下界の幅）✓

定理 4.2 の 2 の仮定（$\|w - P_hw\|_a \le \eta\|Aw\|$、$w \in D(A)$）のもとで、$k \le \dim V_h$ について

$$\frac{\lambda_k^h}{1 + \eta^2\lambda_k^h} \le \lambda_k \le \lambda_k^h$$

が成り立つ。区間の幅は $\dfrac{\eta^2(\lambda_k^h)^2}{1 + \eta^2\lambda_k^h} \le \eta^2(\lambda_k^h)^2$ である。また $\eta^2\lambda_k \lt 1$ なら $\lambda_k^h \le \lambda_k/(1 - \eta^2\lambda_k)$ である。

証明（命題 6.1）

上界は 定理 4.1 です。補題 5.1 から、すべての $w \in V$ について $\|w - P_hw\| \le \eta\|w - P_hw\|_a$ なので、定理 6.3［無限次元の固有値の包含］ の仮定が $C_h = \eta$ で満たされ、$\lambda_k \ge \lambda_k^h/(1 + \eta^2\lambda_k^h)$ です。幅は $\lambda_k^h - \lambda_k^h/(1 + \eta^2\lambda_k^h) = \eta^2(\lambda_k^h)^2/(1 + \eta^2\lambda_k^h)$。下界の式を $\lambda_k^h$ について解くと、$\lambda_k^h(1 - \eta^2\lambda_k) \le \lambda_k$ です。

∎

挟む区間の幅は $O(\eta^2\lambda_k^2)$ で、定理 4.2 の誤差の次数と同じです。つまり、Liu の下界は「大きさの順序が合っているだけの粗い下界」ではなく、$h \to 0$ で Ritz 値と同じ速さで真の値に近づきます。最後の式は 定理 4.2 の 2 の評価の定数を 2 倍改善しています。1 次元の区分 1 次要素（$\eta = h/\pi$）では、幅は $(\lambda_k^h)^2h^2/\pi^2$ 以下で、正確な誤差 $\lambda_k^2h^2/12$（定理 4.3）の約 $1.2$ 倍です。

この挟み込みを証明として使うには、定数 $\eta$（2 次元では補間誤差の定数）の具体的な値と、行列の固有値 $\lambda_k^h$ の計算の丸め誤差の評価が要ります。それが固有値の精度保証の章の内容です。そこでは、より鋭い下界を与える Lehmann–Goerisch の方法（定理 5.2［無限次元の固有値の包含］）も扱っています。

7. 領域の切断

全直線の Schrödinger 作用素 $-\frac{d^2}{dx^2} + V$ の固有値を計算するには、区間 $(-R, R)$ に切って両端で 0 とします。これも Galerkin 近似の一種で（$H^1_0(-R,R)$ の元を 0 で延長すれば全直線の関数になる）、誤差の向きと大きさが同じ道具で分かります。

この節では、$V : \mathbb{R} \to \mathbb{R}$ を連続で下に有界な関数とし、

$$a(u, v) = \int_{\mathbb{R}}(u'v' + Vuv)\,dx,\qquad R(w) = \frac{a(w,w)}{\|w\|^2},\qquad E_k(R) = \inf_{\substack{F \subset H^1_0(-R,R)\\ \dim F = k}}\ \max_{w\in F,\ w\ne0}R(w)$$

とおきます（$\|\cdot\|$ は $L^2(\mathbb{R})$ のノルム）。$E_k(R)$ は、区間 $(-R,R)$ の Dirichlet 問題 $-u'' + Vu = Eu$、$u(\pm R) = 0$ の小さい方から $k$ 番目の固有値です（$V$ に定数を足して正にすれば §2 の設定に入り、補題 2.1 が使えます。定数を足すと $R(w)$ と固有値は同じだけずれます）。

次の定理の 2 は、全直線の固有関数の裾を切り落とし関数で滑らかに 0 に落として箱の中の試験関数にし、その Rayleigh 商が $E_k$ からどれだけ増えるかを測ります。増える分は、固有関数の裾の部分の大きさ $\delta$ だけで押さえられます。証明の中の長い積分の展開は、「エネルギーの増加分がちょうど切り落としの傾きのある薄い層での $u^2$ の積分になる」という等式を出すための計算です。

定理 7.1（領域の切断による固有値の誤差）✓

上の設定で次が成り立つ。

1. $0 \lt R \lt R'$ なら $E_k(R') \le E_k(R)$。
2. 実数 $E_1 \le \cdots \le E_k$ と実数値の $u_1, \ldots, u_k \in H^1(\mathbb{R})$ が、$\langle u_i, u_j\rangle = \delta_{ij}$ と $$\int_{\mathbb{R}}(u_i'v' + Vu_iv)\,dx = E_i\int_{\mathbb{R}}u_iv\,dx\qquad(\text{台がコンパクトなすべての}\ v \in H^1(\mathbb{R}))$$ を満たすとする（全直線の問題の固有値と固有関数）。$E^{*} = \max_i|E_i|$、$R \gt 1$、$\delta = \sum_{i=1}^{k}\int_{|x| \gt R-1}u_i^2\,dx$ とおく。$\delta \lt 1$ ならば $$E_k(R) \le E_k + \frac{(2E^{*} + 1)\,\delta}{1 - \delta}.$$
3. さらに $|u_i(x)| \le Ce^{-\kappa|x|}$（$i \le k$）なら $\delta \le \dfrac{kC^2}{\kappa}e^{-2\kappa(R-1)}$ で、誤差の上界は $O(e^{-2\kappa R})$ である。
4. $E_1, \ldots, E_k$ が、$H^1_0(-R,R)$（0 で延長）を含む空間 $V$ の上の形式 $a$ について 補題 2.1 の $\lambda_1, \ldots, \lambda_k$ であるとき（$V(x) \to \infty$ の場合）、$E_k \le E_k(R)$ である。

証明（定理 7.1）

1. $H^1_0(-R,R) \subset H^1_0(-R',R')$（0 で延長）なので、下限をとる部分空間の範囲が広がります。

2. $\chi$ を、$|x| \le R - 1$ で 1、$R - 1 \le |x| \le R$ で $R - |x|$、$|x| \ge R$ で 0 の関数とします。$\chi$ は Lipschitz 連続で $0 \le \chi \le 1$、$|\chi'| \le 1$、$\chi'$ は $R - 1 \le |x| \le R$ の外で 0 です。$w_i = \chi u_i$ は $H^1_0(-R,R)$ に属し、$w_i' = \chi'u_i + \chi u_i'$ です。

等式。 $a(w_i, w_j) = \frac{E_i + E_j}{2}\langle w_i, w_j\rangle + \int\chi'^2u_iu_j\,dx$ を示します。展開すると

$$a(w_i, w_j) = \int\bigl(\chi'^2u_iu_j + \chi\chi'(u_iu_j)' + \chi^2u_i'u_j' + V\chi^2u_iu_j\bigr)dx$$

です。$u_i$ の式で $v = \chi^2u_j$（台がコンパクトで $H^1$）とすると、$v' = 2\chi\chi'u_j + \chi^2u_j'$ から

$$\int\bigl(2\chi\chi'u_i'u_j + \chi^2u_i'u_j' + V\chi^2u_iu_j\bigr)dx = E_i\langle w_i, w_j\rangle$$

で、$i$ と $j$ を入れ替えた式との平均をとると、$u_i'u_j + u_iu_j' = (u_iu_j)'$ から

$$\int\bigl(\chi\chi'(u_iu_j)' + \chi^2u_i'u_j' + V\chi^2u_iu_j\bigr)dx = \frac{E_i + E_j}{2}\langle w_i, w_j\rangle$$

です。これを展開の式に入れれば等式が出ます。

試験空間。 $c \in \mathbb{R}^k$、$|c| = 1$ について $u = \sum_ic_iu_i$、$\tilde u = \sum_ic_iE_iu_i$、$w = \sum_ic_iw_i = \chi u$ とおきます。等式から

$$a(w, w) = \sum_{i,j}c_ic_j\frac{E_i + E_j}{2}\langle\chi u_i, \chi u_j\rangle + \int\chi'^2u^2 = \int\chi^2\tilde u\,u\,dx + \int\chi'^2u^2\,dx$$

です。$\int\tilde uu = \sum_ic_i^2E_i \le E_k$ です。$1 - \chi^2$ と $\chi'^2$ は $|x| \le R - 1$ で 0、それ以外で 1 以下なので、各点での Cauchy–Schwarz の不等式 $u^2 \le \sum_iu_i^2$、$\tilde u^2 \le (E^{*})^2\sum_iu_i^2$ から

$$\int(1-\chi^2)u^2 \le \delta,\qquad \int\chi'^2u^2 \le \delta,\qquad \Bigl|\int(1 - \chi^2)\tilde uu\Bigr| \le \Bigl(\int_{|x| \gt R-1}\tilde u^2\Bigr)^{1/2}\Bigl(\int_{|x| \gt R-1}u^2\Bigr)^{1/2} \le E^{*}\delta$$

です。したがって

$$a(w,w) \le E_k + (E^{*} + 1)\delta =: N_0,\qquad 1 - \delta \le \|w\|^2 = 1 - \int(1-\chi^2)u^2 \le 1$$

です。とくに $w \ne 0$ なので $w_1, \ldots, w_k$ は 1 次独立で、$F = \mathrm{span}\{w_i\}$ は $H^1_0(-R,R)$ の $k$ 次元部分空間です。$a(w,w) \ge 0$ なら $R(w) \le N_0/(1 - \delta) = E_k + (E_k\delta + (E^{*}+1)\delta)/(1-\delta) \le E_k + (2E^{*}+1)\delta/(1-\delta)$、$a(w,w) \lt 0$ なら $\|w\|^2 \le 1$ から $R(w) \le a(w,w) \le N_0 \le E_k + (2E^{*}+1)\delta/(1-\delta)$ です。$E_k(R) \le \max_FR$ から結論が出ます。

3. $\int_{|x| \gt R-1}u_i^2 \le 2C^2\int_{R-1}^\infty e^{-2\kappa x}dx = \frac{C^2}{\kappa}e^{-2\kappa(R-1)}$ を $i$ について足します。

4. 補題 2.1 の最小値は $V$ のすべての $k$ 次元部分空間についてとるので、$H^1_0(-R,R)$ の中に限った下限 $E_k(R)$ 以下です（定理 4.1 と同じ議論）。

∎

1 と 4 は Dirichlet ラプラシアンの領域の単調性（系 3.1［Dirichlet ラプラシアンの固有値と Weyl 則］）と同じ主張です。2 の証明の等式は、固有関数に切り落とし関数を掛けても、エネルギーの増加が $\int\chi'^2u^2$（切り落としの勾配がある層での固有関数の大きさ）だけで済むことを言っています。誤差が固有関数の裾の 2 乗で決まるのは、定理 4.2 と同じく、固有値の誤差が固有関数の誤差の 2 乗の次数だからです。層の幅を 1 にとったのは式を簡単にするためで、幅 $\ell$ の層を使えば $\chi'^2 \le \ell^{-2}$、裾は $|x| \gt R - \ell$ になります（演習 9.3）。

注意 7.2（本質スペクトルがある場合の下界）

4 は、全直線の作用素が離散スペクトルだけをもつ場合（$V \to \infty$）の主張です。$V \to 0$（$|x| \to \infty$）のように本質スペクトル $[0, \infty)$ がある場合も、その下端より下の固有値については、形式についてのミニマックス原理（Reed–Simon の教科書、参考文献 [3] の第 XIII 章）から同じ不等式 $E_k \le E_k(R)$ が成り立ちます。既存の章のミニマックス原理（定理 6.3［作用素の収束と摂動］）は試験空間が作用素の定義域に含まれる場合の形なので、ここでは定理に含めず、下の例で数値を確かめるにとどめます。基底状態が正にとれて単純であること（定理 5.4［Feller 過程と Feynman–Kac］ が仮定として使う性質。この章では証明しません）は、切断した問題でも同じです。

例 7.3（調和振動子と Pöschl–Teller ポテンシャル）

(a) 調和振動子 $V = x^2$。固有値は $E_k = 2k - 1$、基底状態は $u_1 = \pi^{-1/4}e^{-x^2/2}$ で、$\delta = \int_{|x| \gt R-1}u_1^2 = \mathrm{erfc}(R - 1)$ です。$k = 1$、$E^{*} = 1$ なので、定理 7.1 の上界は $3\delta/(1-\delta)$ です。

(b) Pöschl–Teller ポテンシャル $V = -2/\cosh^2x$。$u_1 = \frac{1}{\sqrt2\cosh x}$ は $-u_1'' + Vu_1 = -u_1$ を満たし（$(1/\cosh x)'' = 1/\cosh x - 2/\cosh^3x$）、$\|u_1\| = 1$、$E_1 = -1$ です。$|u_1| \le \sqrt2e^{-|x|}$（$\kappa = 1$）で、$\delta = \int_{|x| \gt R-1}\frac{dx}{2\cosh^2x} = 1 - \tanh(R-1)$、上界は同じく $3\delta/(1-\delta)$ です。

区間の Dirichlet 固有値は、Chebyshev 選点法（120 点）で求めます。離散化の誤差は表の値より十分小さいです。

# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1
import numpy as np
from scipy.special import erfc

def cheb(N):
    j = np.arange(N + 1)
    x = np.cos(np.pi * j / N)
    c = np.where((j == 0) | (j == N), 2.0, 1.0) * (-1.0) ** j
    D = (c[:, None] / c[None, :]) / (x[:, None] - x[None, :] + np.eye(N + 1))
    D -= np.diag(D.sum(axis=1))
    return D, x

def box_eigs(V, R, N=120, k=3):
    """-u'' + V u = E u（u(±R) = 0）の小さい方から k 個の固有値（Chebyshev 選点法）。"""
    D, x = cheb(N)
    H = -(D @ D)[1:N, 1:N] / R**2 + np.diag(V(R * x[1:N]))
    return np.sort(np.linalg.eigvals(H).real)[:k]

print("harmonic oscillator V = x^2 (E_k = 2k - 1)")
print("    R    E_1(R) - 1    E_2(R) - 3    E_3(R) - 5    bound for k = 1: 3 delta / (1 - delta)")
for R in [2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0]:
    E = box_eigs(lambda x: x**2, R)
    delta = erfc(R - 1)
    print(f"  {R:3.1f}   {E[0] - 1:.3e}     {E[1] - 3:.3e}     {E[2] - 5:.3e}     {3 * delta / (1 - delta):.3e}")

print("Poschl-Teller V = -2 / cosh^2 x (E_1 = -1)")
print("    R    E_1(R) + 1    bound 3 delta / (1 - delta)    ratio of successive errors (R -> R + 1)")
prev = None
for R in [3.0, 4.0, 5.0, 6.0, 7.0, 8.0]:
    E = box_eigs(lambda x: -2 / np.cosh(x) ** 2, R, k=1)[0] + 1
    delta = 1 - np.tanh(R - 1)
    print(f"  {R:3.1f}   {E:.3e}     {3 * delta / (1 - delta):.3e}" + ("" if prev is None else f"                      {prev / E:.3f}"))
    prev = E
print("  e^2 =", f"{np.exp(2):.3f}")

harmonic oscillator V = x^2 (E_k = 2k - 1)
    R    E_1(R) - 1    E_2(R) - 3    E_3(R) - 5    bound for k = 1: 3 delta / (1 - delta)
  2.0   7.492e-02     5.296e-01     1.800e+00     5.600e-01
  2.5   9.908e-03     1.028e-01     4.734e-01     1.053e-01
  3.0   7.822e-04     1.216e-02     8.225e-02     1.410e-02
  3.5   3.609e-05     7.990e-04     7.940e-03     1.221e-03
  4.0   9.817e-07     2.921e-05     4.024e-04     6.627e-05
  4.5   1.588e-08     6.083e-07     1.101e-05     2.229e-06
  5.0   1.535e-10     7.343e-09     1.680e-07     4.625e-08
Poschl-Teller V = -2 / cosh^2 x (E_1 = -1)
    R    E_1(R) + 1    bound 3 delta / (1 - delta)    ratio of successive errors (R -> R + 1)
  3.0   2.025e-02     1.119e-01
  4.0   2.695e-03     1.491e-02                      7.514
  5.0   3.635e-04     2.013e-03                      7.414
  6.0   4.916e-05     2.724e-04                      7.394
  7.0   6.652e-06     3.687e-05                      7.390
  8.0   9.003e-07     4.989e-06                      7.389
  e^2 = 7.389

誤差はすべて正で、$R$ とともに単調に減り、上界の内側にあります。Pöschl–Teller では $R$ を 1 増やすごとに誤差が $1/e^2$ になり、定理 7.1 の 3 の $e^{-2\kappa R}$（$\kappa = 1$）のとおりです。調和振動子では固有関数が Gauss 関数的に減衰するので、誤差は $e^{-R^2}$ のように、指数関数より速く減ります。励起状態ほど固有関数が外に広がるので、同じ $R$ での誤差は大きくなります。必要な箱の大きさは、求めたい最も高い固有状態の広がりで決まります。

8. 偽の固有値とスペクトル法

8.1. 偽の固有値が現れない範囲

離散化した行列の固有値が、作用素のスペクトルと無関係な場所に現れる現象をスペクトル汚染、その固有値を偽の固有値といいます。§2 の設定（離散スペクトルだけをもつ正の作用素）では、定理 4.1 から $k$ 番目の Ritz 値は $k$ 番目の固有値より下に来ないので、下からの汚染はありません。本質スペクトルをもつ作用素でも、その下端より下では同じです。

命題 8.1（本質スペクトルの下端より下には偽の固有値が現れない）✓

$\mathcal{H}$ を無限次元の Hilbert 空間、$A$ を下に有界な自己共役作用素、$\Sigma = \inf\sigma_{\mathrm{ess}}(A)$、$\mu_k(A)$ をミニマックス値（定義 6.1［作用素の収束と摂動］）とする。$V_h \subset D(A)$ を $n$ 次元の部分空間、$\Lambda_1 \le \cdots \le \Lambda_n$ を、$V_h$ の正規直交基底 $\psi_i$ についての行列 $(\langle\psi_i, A\psi_j\rangle)$ の固有値（Ritz 値）とする。

1. $k = 1, \ldots, n$ について $\Lambda_k \ge \mu_k(A)$。
2. $c \lt \Sigma$ とすると、$c$ より小さい Ritz 値の個数は、$c$ より小さい $A$ の固有値の個数（重複度こみ）以下である。とくに $\inf\sigma(A)$ より小さい Ritz 値はなく、$\Lambda_k \lt \Sigma$ ならば $A$ は $\Sigma$ より小さい固有値を $k$ 個以上もち、小さい方から $k$ 番目の固有値は $\Lambda_k$ 以下である。

証明（命題 8.1）

1. 行列 $(\langle\psi_i, A\psi_j\rangle)$ はエルミートで、Courant–Fischer の定理（定理 3.1［行列の固有値評価と不等式］。固有値を小さい順に数え直す）から、$V_h$ の $k$ 次元部分空間 $F$ で $\max\{\langle\psi, A\psi\rangle : \psi \in F,\ \|\psi\| = 1\} = \Lambda_k$ となるものがあります。$F \subset D(A)$ なので、系 6.4［作用素の収束と摂動］ から $\mu_k(A) \le \Lambda_k$ です。

2. $c$ より小さい Ritz 値が $k$ 個あるとすると $\Lambda_k \lt c \lt \Sigma$ で、1 から $\mu_k(A) \lt \Sigma$ です。ミニマックス原理（定理 6.3［作用素の収束と摂動］）の 2 つの場合のうち、$\mu_k = \Sigma$ となる場合は起こらないので、$A$ は $\Sigma$ より小さい固有値を $k$ 個以上もち、$\mu_k(A)$ は小さい方から $k$ 番目の固有値です。$\mu_1 \le \cdots \le \mu_k \le \Lambda_k \lt c$ なので、$c$ より小さい固有値が $k$ 個以上あります。$\mu_1(A) = \inf\sigma(A) \le \Lambda_1$ から残りの主張が出ます。

∎

次の例の直感は単純です。この作用素は、一方の成分には $+1$ 倍、もう一方には $-1$ 倍として働きます。両方の成分を混ぜた試験ベクトルの Rayleigh 商は、$+1$ と $-1$ の重みつき平均になるので、そのあいだの値がいくらでも現れます。

例 8.2（スペクトルのすき間に現れる偽の固有値）

$\mathcal{H} = \ell^2(\mathbb{N})\oplus\ell^2(\mathbb{N})$、$A(x \oplus y) = x \oplus (-y)$ とします。$A$ は有界な自己共役作用素で、スペクトルは $\{-1, 1\}$（どちらも無限重の固有値）、$\Sigma = -1$ です。$e$ を $\ell^2(\mathbb{N})$ の単位ベクトルとし、$\psi_t = \cos t\,(e \oplus 0) + \sin t\,(0 \oplus e)$ の張る 1 次元の部分空間をとると、Ritz 値は $\langle\psi_t, A\psi_t\rangle = \cos^2t - \sin^2t = \cos2t$ です。$t$ を動かすと $[-1, 1]$ のどの値にもなり、すき間 $(-1, 1)$ の中の値は $A$ のスペクトルに属しません。命題 8.1 が保証するのは $\Sigma = -1$ より下だけで、本質スペクトルの 2 つの部分に挟まれたすき間では、試験関数が両側の成分を混ぜると偽の固有値が現れます。Dirac 作用素（正と負のエネルギーの連続スペクトルのすき間に束縛状態がある）や、周期的な媒質の欠陥の準位の計算で実際に起こります。

注意 8.3（すき間の中の固有値（文献））

すき間の中で偽の固有値が現れる条件と、それを避ける方法（試験空間の選び方、2 次の相対スペクトルなど）は、Davies–Plum と Lewin–Séré の論文（参考文献 [6]、[7]）にあります。この章では扱いません。

8.2. スペクトル法の固有値

定理 4.2 の 1 は、近似空間が何であっても成り立ちます。近似空間を三角多項式にとれば、固有関数の滑らかさに応じた速さで収束します。

命題 8.4（Fourier–Galerkin 法の固有値の指数収束）✓

$\mathcal{H}$ を周期 $2\pi$ の実数値関数の $L^2$ 空間（内積 $\langle u, v\rangle = \frac{1}{2\pi}\int_{-\pi}^{\pi}uv\,dx$）、$V_c$ を連続な周期関数で $1 \le V_c(x) \le V_1$ を満たすものとし、$V = H^1_{\mathrm{per}}$ 上の形式 $a(u,v) = \langle u', v'\rangle + \langle V_cu, v\rangle$ が §2 の設定を満たすとする（固有値 $\lambda_m$、固有関数 $u_m$）。$N = 2K+1$、$V_h = S_N$（次数 $K$ 以下の実三角多項式の空間。スペクトル法の 定義 2.1［スペクトル法］）とし、その Ritz 値を $\lambda_k^N$ と書く。$u_1, \ldots, u_k$ が帯状領域 $\{|\mathrm{Im}\,z| \lt a_0\}$ に正則に延長され、そこで $|u_i| \le M_0$ を満たすなら、$0 \lt b \lt a_0$ について、$N$ によらない定数 $C$ があって、十分大きいすべての $N$ で

$$0 \le \lambda_k^N - \lambda_k \le C\,(N+1)^2e^{-bN}$$

が成り立つ。$C = 4kM_0^2G_be^{b}V_1\lambda_k/\lambda_1$、$G_b = \sum_{j\ge1}j^2e^{-2bj}$ ととれる。

証明（命題 8.4）

$u = \sum_{i\le k}c_iu_i$、$\sum c_i^2 = 1$ とします。Ritz 射影は $a$ のノルムでの最良近似なので、Fourier 切断 $P_Nu \in S_N$ と比べて $\|u - P_hu\|_a \le \|u - P_Nu\|_a$ です。$a(v,v) \le V_1(\|v'\|^2 + \|v\|^2) = V_1\|v\|_{H^1}^2$ です。Fourier 係数は 補題 3.1［数値積分 II］ と Cauchy–Schwarz の不等式から $|\hat u(n)| \le \sum_i|c_i||\hat u_i(n)| \le \sqrt kM_0e^{-b|n|}$ なので

$$\|u - P_Nu\|_{H^1}^2 = \sum_{|n| \gt K}(1 + n^2)|\hat u(n)|^2 \le 4kM_0^2\sum_{n \gt K}n^2e^{-2bn} \le 4kM_0^2(K+1)^2e^{-2bK}\sum_{j\ge1}j^2e^{-2bj}$$

です（$1 + n^2 \le 2n^2$ と、$n = K + j \le (K+1)j$ を使った）。よって 定理 4.2 の $\varepsilon_k$ は $\varepsilon_k^2 \le 4kM_0^2G_bV_1(K+1)^2e^{-2bK}$ を満たし、$N \to \infty$ で 0 に収束します。$4\varepsilon_k^2 \le \lambda_1$ となる $N$ では、同じ定理の 1 から

$$0 \le \lambda_k^N - \lambda_k \le \lambda_k\cdot\frac{2\varepsilon_k^2/\lambda_1}{1 - 2\varepsilon_k^2/\lambda_1} \le \frac{4\lambda_k}{\lambda_1}\varepsilon_k^2$$

です。$(K+1)^2 = (N+1)^2/4$、$e^{-2bK} = e^{b}e^{-bN}$ を入れると主張の式が出ます。

∎

例 8.5（Mathieu 型の問題の固有値）

$-u'' + 2q\cos(2x)\,u = \lambda u$（周期 $2\pi$）の固有値は、Mathieu 関数の特性値 $a_m(q)$、$b_m(q)$ です。ポテンシャルに定数を足せば 命題 8.4 の形になり（固有値の差は変わらない）、係数が整関数の線形常微分方程式の解は整関数なので、固有関数は任意の幅の帯状領域で正則です。$e^{inx}$ を基底にとると、Galerkin 法の行列は $\mathrm{diag}(n^2)$ に、$n$ と $n \pm 2$ を結ぶ成分 $q$ を加えたものです。$q = 5$ で、下位 3 個の固有値を Fourier–Galerkin 法と Fourier 選点法（スペクトル法の 2 階のスペクトル微分行列に $\mathrm{diag}(2q\cos2x_j)$ を加えた行列）で求め、scipy.special の特性値と比べます。

# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1
import numpy as np
from scipy.special import mathieu_a, mathieu_b

q = 5.0
ref = np.sort([mathieu_a(0, q), mathieu_b(1, q), mathieu_a(1, q)])
print("reference (scipy.special):", "  ".join(f"{r:.10f}" for r in ref))
print("    N    Galerkin: lambda_k^N - lambda_k (k = 1, 2, 3)      collocation: lambda_k^N - lambda_k")
for N in [9, 13, 17, 21, 25, 29]:
    K = (N - 1) // 2
    n = np.arange(-K, K + 1)
    G = np.diag(n.astype(float) ** 2) + q * (np.eye(N, k=2) + np.eye(N, k=-2))   # e^{inx} についての行列
    gal = np.linalg.eigvalsh(G)[:3]
    x = 2 * np.pi * np.arange(N) / N
    k = np.fft.fftfreq(N, d=1.0 / N)
    D2 = np.real(np.fft.ifft(-(k**2)[:, None] * np.fft.fft(np.eye(N), axis=0), axis=0))
    col = np.linalg.eigvalsh(-D2 + np.diag(2 * q * np.cos(2 * x)))[:3]
    print(f"{N:5d}    " + "  ".join(f"{l - r:+.2e}" for l, r in zip(gal, ref)) + "        "
          + "  ".join(f"{l - r:+.2e}" for l, r in zip(col, ref)))

reference (scipy.special): -5.8000460209  -5.7900805986  1.8581875415
    N    Galerkin: lambda_k^N - lambda_k (k = 1, 2, 3)      collocation: lambda_k^N - lambda_k
    9    +1.28e-02  +8.95e-02  +3.32e-01        -4.50e-02  -4.18e-02  -3.04e-01
   13    +1.13e-04  +1.39e-03  +1.02e-02        -1.64e-04  -4.91e-04  -4.75e-03
   17    +3.88e-07  +7.44e-06  +7.00e-05        -3.75e-07  -1.89e-06  -2.02e-05
   21    +6.15e-10  +1.70e-08  +1.87e-07        -4.52e-10  -3.43e-09  -4.08e-08
   25    +4.88e-13  +1.92e-11  +2.36e-10        -3.12e-13  -3.26e-12  -4.18e-11
   29    +8.88e-16  +9.77e-15  +1.44e-13        +2.22e-14  +0.00e+00  -8.84e-14

$N$ を 4 増やすごとに誤差が 2〜3 桁減り、$N = 29$ で丸めの大きさです。2 次の中心差分で同じ精度（$10^{-13}$）を得るには、第 1 固有値の誤差が約 $1.0\,h^2$ なので（$N = 200, 400, 800$ 点の周期的な 3 点差分で確かめた値）、$h \approx 3\times10^{-7}$、すなわち $10^7$ 個程度の格子点が要ることになります。Galerkin 法の誤差はすべて正です（定理 4.1）。選点法の誤差は同じ速さで減りますが、符号は負で、上界にはなっていません。選点法は積 $V_cu$ を格子点で標本化するので折り返しが入り（定理 6.3［スペクトル法］）、適合な Galerkin 近似ではないからです。

9. 演習

演習 9.1易

区間のモデル問題を $h = 1/4$ で離散化する。3 点差分の第 1 固有値、区分 1 次要素の第 1 Ritz 値、両者の平均を求め、$\lambda_1 = \pi^2 \approx 9.8696$ と比べよ。

解答

$\theta_1 = \pi/4$ です。定理 3.1 から $\lambda_1^{\mathrm{FD}} = 64\sin^2(\pi/8) = 32(1 - \cos(\pi/4)) \approx 9.3726$。定理 4.3 から $\lambda_1^h = 96\,\frac{1 - \cos(\pi/4)}{2 + \cos(\pi/4)} \approx 96\times\frac{0.292893}{2.707107} \approx 10.3866$。誤差は $-0.497$ と $+0.517$ で、先頭項 $\lambda_1^2h^2/12 \approx 0.507$ に近い値です。平均は $9.8796$ で、誤差は $+0.0100$。予言値 $\lambda_1^3h^4/360 \approx 961.4/(256\times360) \approx 0.0104$ と合います。未知数 3 個の行列 2 つから、相対誤差 $0.1\%$ の値が得られています。

演習 9.2標準

§2 の設定で、$\lambda_k$ の固有ベクトル $u_k$（$\|u_k\| = 1$）と、$\|v\| = 1$ である任意の $v \in V$ について

$$R(v) - \lambda_k = \|v - u_k\|_a^2 - \lambda_k\|v - u_k\|^2$$

を示せ。これから、$R(v)$ の $\lambda_k$ からのずれが $v$ の誤差の 2 乗の次数であることを説明せよ。

解答

$a(u_k, v) = \lambda_k\langle u_k, v\rangle$、$a(u_k,u_k) = \lambda_k$ から $\|v - u_k\|_a^2 = a(v,v) - 2\lambda_k\langle u_k, v\rangle + \lambda_k$ です。$\|v - u_k\|^2 = 2 - 2\langle u_k, v\rangle$ なので $\lambda_k\|v - u_k\|^2 = 2\lambda_k - 2\lambda_k\langle u_k, v\rangle$。差をとると $a(v,v) - \lambda_k = R(v) - \lambda_k$ です。右辺は $\|v - u_k\|_a^2$ 以下なので、$v$ がエネルギーノルムで $u_k$ に $\varepsilon$ まで近ければ、Rayleigh 商は $\lambda_k$ に $\varepsilon^2$ まで近くなります。命題 5.2 のエネルギーノルムの等式は、$v = u_k^h$ の場合です。$k = 1$ なら $\|w\|_a^2 \ge \lambda_1\|w\|^2$ から右辺は 0 以上で、Rayleigh–Ritz の上界 $R(v) \ge \lambda_1$ が再び出ます。$k \ge 2$ では右辺は負にもなり、1 つの試験関数の Rayleigh 商は $\lambda_k$ の上界とは限りません。

演習 9.3標準

定理 7.1 の設定で $k = 1$ とし、$0 \lt \ell \lt R$、$\delta_\ell = \int_{|x| \gt R - \ell}u_1^2\,dx \lt 1$ とする。幅 $\ell$ の層で 1 から 0 に 1 次式で落ちる切り落とし関数を使って

$$E_1(R) - E_1 \le \frac{\delta_\ell}{\ell^2(1 - \delta_\ell)}$$

を示せ。調和振動子の基底状態で $R = 4$、$\ell = 1$ のときの値を、例 7.3 の表と比べよ。

解答

$\chi$ を、$|x| \le R - \ell$ で 1、$R - \ell \le |x| \le R$ で $(R - |x|)/\ell$、それ以外で 0 とします。定理 7.1 の証明の等式は $\chi$ の形によらず成り立ち、$i = j = 1$ で $a(\chi u_1, \chi u_1) = E_1\|\chi u_1\|^2 + \int\chi'^2u_1^2$ です。$\chi'^2 \le \ell^{-2}$ で、$\chi'$ は $|x| \gt R - \ell$ の外で 0 なので $\int\chi'^2u_1^2 \le \delta_\ell/\ell^2$、また $\|\chi u_1\|^2 \ge 1 - \delta_\ell$ です。$\chi u_1 \in H^1_0(-R,R)$ の張る 1 次元の空間を試験空間にとると

$$E_1(R) \le R(\chi u_1) = E_1 + \frac{\int\chi'^2u_1^2}{\|\chi u_1\|^2} \le E_1 + \frac{\delta_\ell}{\ell^2(1 - \delta_\ell)}$$

です。$k = 1$ では等式が正確なので、定理の一般の評価にある因子 $2E^{*} + 1$ が要りません。調和振動子で $R = 4$、$\ell = 1$ なら $\delta_1 = \mathrm{erfc}(3) \approx 2.21\times10^{-5}$ で、上界は $2.21\times10^{-5}$。表の実際の誤差は $9.8\times10^{-7}$、定理の上界は $6.6\times10^{-5}$ でした。$\ell$ を大きくすると $\ell^{-2}$ は減り $\delta_\ell$ は増えるので、最良の $\ell$ はその釣り合いで決まります。

演習 9.4難

1 次元の等間隔の区分 1 次要素で、質量行列を集中化した固有値問題 $K\boldsymbol{x} = \lambda hI\boldsymbol{x}$ は、内積を $\langle u, v\rangle_h = h\sum_ju(x_j)v(x_j)$（台形則）に取り替えた Galerkin 近似とみなせることを示せ。$v \in V_h$ について $\langle v, v\rangle_h \ge \|v\|^2$ であることを示し、これから集中質量の固有値が一致質量の Ritz 値以下になること（$\lambda_k^{\mathrm{FD}} \le \lambda_k^h$）を、定理 3.1 と 定理 4.3 の閉じた形を使わずに導け。

解答

$u = \sum_jx_j\varphi_j \in V_h$ について $u(x_j) = x_j$ なので、$\langle\varphi_j, \varphi_i\rangle_h = h\delta_{ij}$ で、$a(u, v) = \lambda\langle u, v\rangle_h$（$v \in V_h$）は $K\boldsymbol{x} = \lambda hI\boldsymbol{x}$ です。

1 つの要素 $[x_j, x_{j+1}]$ の上で $v$ は 1 次式で、端の値を $\alpha, \beta$ とすると $\int v^2 = \frac h3(\alpha^2 + \alpha\beta + \beta^2)$、台形則は $\frac h2(\alpha^2 + \beta^2)$ です。差は $\frac h6(\alpha^2 - 2\alpha\beta + \beta^2) = \frac h6(\alpha - \beta)^2 \ge 0$。要素について足すと $\langle v, v\rangle_h \ge \|v\|^2$ です（$v(0) = v(1) = 0$ なので、台形則の和は $h\sum_jv(x_j)^2$ に等しい）。

補題 2.3 の証明は内積を $\langle\cdot,\cdot\rangle_h$ に替えても通るので、集中質量の第 $k$ 固有値は $\min_{\dim F = k}\max_{v\in F}a(v,v)/\langle v, v\rangle_h$ です。各 $v$ で $a(v,v)/\langle v,v\rangle_h \le a(v,v)/\|v\|^2$ なので、$\max$ も $\min$ も大きくならず、$\lambda_k^{\mathrm{FD}} \le \lambda_k^h$ です。差は $\frac h6\sum_j(v(x_{j+1}) - v(x_j))^2 = \frac{h^2}{6}\|v'\|^2$ で、$\langle v, v\rangle_h = \|v\|^2 + \frac{h^2}{6}a(v,v)$ という等式になっています。これを Rayleigh 商に入れると $\frac{1}{\lambda^{\mathrm{FD}}} = \frac{1}{\lambda^h} + \frac{h^2}{6}$ が同じ固有ベクトルについて成り立ち、定理 4.3 の 2 つの閉じた形の関係（$\frac{h^2}{4\sin^2(\theta/2)} - \frac{h^2(2 + \cos\theta)}{6(1 - \cos\theta)} = \frac{h^2}{6}$）と一致します。集中化は、内積を大きくすることで固有値を下げます。真の値 $\lambda_k$ を下回るかどうかはこの議論からは分からず、定数係数の場合は 定理 3.1 の計算で確かめられた事実です。

参考文献

1. G. Strang, G. J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, 1973 — 第 6 章（固有値問題。Ritz 値の誤差が固有関数の誤差の 2 乗になること）、第 8 章（特異点のある領域）。
2. I. Babuška, J. Osborn, “Eigenvalue problems”, in Handbook of Numerical Analysis, Vol. II (P. G. Ciarlet, J. L. Lions, eds.), North-Holland, 1991, 641–787.
3. M. Reed, B. Simon, Methods of Modern Mathematical Physics IV: Analysis of Operators, Academic Press, 1978 — 第 XIII 章（ミニマックス原理）。
4. D. Boffi, “Finite element approximation of eigenvalue problems”, Acta Numerica 19 (2010), 1–120.
5. T. Betcke, L. N. Trefethen, “Reviving the method of particular solutions”, SIAM Review 47 (2005), 469–491 — L 字形領域の固有値の高精度の値。
6. E. B. Davies, M. Plum, “Spectral pollution”, IMA Journal of Numerical Analysis 24 (2004), 417–438.
7. M. Lewin, É. Séré, “Spectral pollution and how to avoid it”, Proceedings of the London Mathematical Society 100 (2010), 864–900.