スペクトル法：Fourier・Chebyshev 選点法と指数収束

0. この記事の要点

関数の値を格子点でしか知らないとき、その微分をどう求めればよいでしょうか。差分法は隣の数点だけを使いますが、この記事で扱うスペクトル法は、すべての格子点の値を使って、波の形をした関数（sin や cos、または多項式）の和で関数全体を近似します。関数が滑らかなときは、点を少し増やすだけで誤差が急に小さくなります。読み終えると、その理由と、うまくいかない場合（不連続な関数、折り返し）の見分け方が分かり、差分法との使い分けを判断できます。

• スペクトル法は、解を区間全体で定義された滑らかな関数（三角多項式・代数多項式）の和で近似します。Fourier 切断 $P_Nu$ の誤差は捨てた Fourier 係数だけで決まり、$u \in H^s$ なら $O(N^{-s})$、$u$ が実軸のまわりの帯状領域で正則なら $O(e^{-bN/2})$ です。差分法や有限要素法のように「次数 $p$ の方法」という上限がありません。
• 等間隔の $N$ 点で関数を標本化すると、波数 $n$ と $n + mN$ の指数関数は区別できません（折り返し）。離散 Fourier 係数は $\tilde u_k = \sum_{m}\hat u(k + mN)$ です。三角補間 $I_Nu$ の誤差は、切断の誤差に折り返しの分が加わったもので、次数は切断と同じです。
• 補間した三角多項式を微分して格子点で値をとる操作は行列 $D$（スペクトル微分行列）で書けます。$D$ は歪対称な巡回行列で、固有値は $ik$（$|k| \le K$）です。2 階微分のスペクトル半径（固有値の絶対値の最大）は $K^2 \approx N^2/4$ なので、陽的な時間積分の刻み幅は $\Delta t = O(N^{-2})$ に制限されます。
• Galerkin 法は残差を近似空間に直交させ、選点法は格子点で方程式を満たさせます。両者の解の差は、右辺と「係数 × 解」の折り返し誤差で抑えられます。2 次の非線形項の折り返しは、$3K + 1$ 点以上の細かい格子で積を作れば消せます（3/2 則）。
• 周期的でない問題では Chebyshev–Gauss–Lobatto 点 $x_j = \cos(j\pi/N)$ を使います。Bernstein 楕円 $E_\rho$（焦点が $\pm1$ で区間 $[-1,1]$ を内側に含む楕円。この内側で関数が正則（複素微分できる）なら誤差が指数的に減る）で正則な関数の補間誤差は $O(\rho^{-N})$、導関数の誤差は $O(N^2\rho^{-N})$ です。点が端に $O(N^{-2})$ の間隔で集まるため、微分行列のノルムは $N^2$、2 階では $N^4$ の大きさになります。
• 不連続な関数では、切断の誤差は一様には 0 に収束せず、不連続点から離れた点でも $O(1/N)$ にとどまります。スペクトル法は滑らかな解に向く方法です。

この記事で使う事前知識

• Fourier 係数: 周期関数を $e^{ikx}$ の和で表したときの、各波の大きさを表す数。→ Fourier 係数
• Parseval の等式: 関数の大きさ（$L^2$ ノルム）の 2 乗が、Fourier 係数の絶対値の 2 乗の和に等しいという関係。→ フーリエ級数
• 正規直交基底: 長さ 1 で互いに直交するベクトル（関数）の集まりで、すべての元がその（無限和を含む）和で書けるもの。→ 正規直交基底
• 多項式補間と Chebyshev 点: 与えられた点を通る多項式を作る方法。等間隔の点では Runge 現象で不安定になり、端に密な Chebyshev 点が安定。→ 多項式補間
• Galerkin 近似: 方程式の残りを近似空間に直交させて近似解を決める方法。→ 有限要素法 I
• スペクトル半径: 行列の固有値の絶対値の最大。→ スペクトル半径

1. 動機

2 次の中心差分で $u'(x)$ を近似すると誤差は $O(h^2)$、4 次なら $O(h^4)$ です。次数を上げるには、より多くの格子点の値を使います。では、すべての格子点の値を使ったらどうなるでしょうか。周期関数の場合、$N$ 個の格子点の値をすべて通る三角多項式を作り、それを微分することになります。これが Fourier 選点法（擬スペクトル法）です。

§4 の数値実験では、$u(x) = e^{\sin x}$ の導関数を $N = 33$ 点で計算すると、2 次の差分の誤差が $2.4 \times 10^{-2}$、4 次の差分が $1.0 \times 10^{-3}$ であるのに対し、スペクトル微分の誤差は $3.2 \times 10^{-15}$ です。誤差が $N$ の冪ではなく指数関数で減るからです。一方 $|\sin x|^3$ のように滑らかさが有限の関数では、3 つの方法の誤差はどれも $N^{-2}$ に比例して減り、差は定数倍にとどまります。

この章の目的は、この振る舞いを定理として述べて証明することです。鍵は 2 つあります。1 つは、近似の誤差が Fourier 係数の減り方で決まること（フーリエ級数の Parseval の等式）。もう 1 つは、格子点での標本化が高い波数を低い波数に重ねてしまうこと（折り返し）で、これは周期関数の台形則が指数的に収束する理由（補外と Romberg 積分）と同じ現象です。

離散 Fourier 変換を $O(N\log N)$ の演算量で計算する算法（FFT）は、この章では道具として使うだけにします。コードでは numpy.fft を呼びます。

2. Fourier 切断と周期 Sobolev 空間

$\mathbb{T} = \mathbb{R}/2\pi\mathbb{Z}$ 上の関数（周期 $2\pi$ の関数）を扱います。記号は フーリエ級数に合わせ、$e_k(x) = e^{ikx}$、

$$\hat u(k) = \frac{1}{2\pi}\int_{-\pi}^{\pi}u(x)e^{-ikx}\,dx,\qquad \langle u, v\rangle = \frac{1}{2\pi}\int_{-\pi}^{\pi}u\,\bar v\,dx,\qquad \|u\| = \langle u, u\rangle^{1/2}$$

とします（定義 2.1［フーリエ級数］）。$\{e_k\}$ は $L^2(\mathbb{T})$ の正規直交基底で、$\|u\|^2 = \sum_k|\hat u(k)|^2$ です（定理 5.1［フーリエ級数］）。$\|u\|_\infty$ は最大値ノルムです。

この章を通じて、格子点の数 $N$ は奇数とし、$N = 2K + 1$ と書きます（偶数の場合の扱いは 注意 5.3 で述べます）。

定義 2.1（周期 Sobolev 空間と Fourier 切断）↗

$s \ge 0$ と $u \in L^2(\mathbb{T})$ について

$$|u|_{H^s}^2 = \sum_{k \ne 0}|k|^{2s}|\hat u(k)|^2,\qquad \|u\|_{H^s}^2 = \sum_{k\in\mathbb{Z}}(1 + k^2)^s|\hat u(k)|^2$$

とおき、$\|u\|_{H^s} \lt \infty$ である $u$ の全体を $H^s_{\mathrm{per}}$ と書く。$S_N = \mathrm{span}\{e_k : |k| \le K\}$（次数 $K$ 以下の三角多項式の空間。次元は $N$）への $L^2$ 直交射影

$$P_Nu = \sum_{|k| \le K}\hat u(k)\,e_k$$

を Fourier 切断という。$u \in H^1_{\mathrm{per}}$ について、Fourier 係数が $ik\hat u(k)$ である $L^2$ の元を $u'$ と書く。

$H^s$ は「$s$ 回微分できて、その導関数の大きさが有限」という滑らかさを、Fourier 係数の減り方で測るものです。$|k|^{2s}$ の重みは、$s$ 回微分すると $k$ 番目の係数が $|k|^s$ 倍になることに対応し、$s$ が整数でなくても同じ式で滑らかさの度合いを表せます。

$u$ が $C^1$ 級なら、部分積分で $\widehat{u'}(k) = ik\hat u(k)$ となるので、この $u'$ はふつうの導関数です。整数 $s$ では $|u|_{H^s} = \|u^{(s)}\|$ です。全空間の Sobolev 空間（定義 3.1［Sobolev 空間］）で Fourier 変換の重み $\langle\xi\rangle^{2s}$ を使ったのと同じ定義を、Fourier 係数について行ったものです。

以下では $s \gt 1/2$ について $Z_s = \sum_{m=1}^{\infty}m^{-2s}$ と書きます。$m \ge 2$ で $m^{-2s} \le \int_{m-1}^{m}x^{-2s}dx$ なので $Z_s \le 1 + \int_1^\infty x^{-2s}dx = 2s/(2s-1)$ です。

補題 2.2（Fourier 係数が絶対総和可能な関数）✓

$u \in L^2(\mathbb{T})$ が $\sum_k|\hat u(k)| \lt \infty$ を満たすなら、$u$ はほとんどいたるところで連続関数 $g(x) = \sum_k\hat u(k)e^{ikx}$（一様収束）に等しい。$s \gt 1/2$ で $u \in H^s_{\mathrm{per}}$ ならこの仮定が成り立ち、

$$\sum_{k \ne 0}|\hat u(k)| \le (2Z_s)^{1/2}\,|u|_{H^s}$$

である。

証明（補題 2.2）

$|\hat u(k)e^{ikx}| = |\hat u(k)|$ で $\sum_k|\hat u(k)| \lt \infty$ なので、級数は一様収束し（系 6.4［関数列と一様収束］）、和 $g$ は連続です（定理 4.1［関数列と一様収束］）。一様収束する級数は項別に積分できるので（定理 5.2［関数列と一様収束］ を実部と虚部に使う）、$\hat g(n) = \sum_k\hat u(k)\langle e_k, e_n\rangle = \hat u(n)$ です。したがって $u - g$ の Fourier 係数はすべて 0 で、Parseval の等式（定理 5.1［フーリエ級数］）から $\|u - g\| = 0$、すなわちほとんどいたるところで $u = g$ です。

後半は Cauchy–Schwarz の不等式（命題 2.2［Hilbert 空間の基礎］ を $\ell^2$ に使う）から

$$\sum_{k\ne0}|\hat u(k)| = \sum_{k\ne0}|k|^{-s}\cdot|k|^{s}|\hat u(k)| \le \Bigl(\sum_{k\ne0}|k|^{-2s}\Bigr)^{1/2}\Bigl(\sum_{k\ne0}|k|^{2s}|\hat u(k)|^2\Bigr)^{1/2} = (2Z_s)^{1/2}|u|_{H^s}$$

です。

∎

以下、この補題の仮定を満たす $u$ は連続な代表元 $g$ と同一視します。全空間での対応する結果は 定理 4.1［Sobolev 空間］ です。

定理 2.3（Fourier 切断の誤差）✓

$0 \le r \le s$、$u \in H^s_{\mathrm{per}}$、$N = 2K + 1$ とする。このとき

$$|u - P_Nu|_{H^r} \le (K+1)^{r-s}\,|u|_{H^s} \le \Bigl(\frac N2\Bigr)^{r-s}|u|_{H^s},\qquad \|u - P_Nu\|_{H^r} \le 2^{r/2}(K+1)^{r-s}\,|u|_{H^s}$$

が成り立つ。$r = 0$ では左辺は $\|u - P_Nu\|$ である。また $P_Nu$ は、$S_N$ の元のうち $\|u - v\|_{H^r}$ を最小にするものである。

証明（定理 2.3）

$u - P_Nu$ の Fourier 係数は、$|k| \le K$ で 0、$|k| \gt K$ で $\hat u(k)$ です。$|k| \ge K + 1$ と $r - s \le 0$ から $|k|^{2r} = |k|^{2(r-s)}|k|^{2s} \le (K+1)^{2(r-s)}|k|^{2s}$ なので

$$|u - P_Nu|_{H^r}^2 = \sum_{|k| \gt K}|k|^{2r}|\hat u(k)|^2 \le (K+1)^{2(r-s)}\sum_{|k| \gt K}|k|^{2s}|\hat u(k)|^2 \le (K+1)^{2(r-s)}|u|_{H^s}^2$$

です。$K + 1 \gt N/2$ から 2 番目の不等式が出ます。$|k| \ge 1$ では $1 + k^2 \le 2k^2$ なので $\|u - P_Nu\|_{H^r}^2 \le 2^r|u - P_Nu|_{H^r}^2$ です。$r = 0$ では $u - P_Nu$ に $k = 0$ の成分がないので $|u - P_Nu|_{H^0} = \|u - P_Nu\|$ です。

最後の主張を示します。$v = \sum_{|k|\le K}b_ke_k \in S_N$ について

$$\|u - v\|_{H^r}^2 = \sum_{|k|\le K}(1+k^2)^r|\hat u(k) - b_k|^2 + \sum_{|k| \gt K}(1+k^2)^r|\hat u(k)|^2$$

で、第 1 項は $b_k = \hat u(k)$ のとき 0 になり、第 2 項は $v$ によりません。

∎

この定理が「スペクトル精度」の中身です。$u$ が滑らかであるほど（$s$ が大きいほど）誤差の次数が上がり、方法の側に次数の上限がありません。$u$ が $C^\infty$ 級なら、誤差は $N$ のどの冪よりも速く減ります。

3. 三角補間と折り返し

切断 $P_Nu$ を作るには Fourier 係数の積分が要ります。計算では、格子点での値だけから近似を作ります。

定義 3.1（離散 Fourier 係数と三角補間）↗

$N = 2K+1$、格子点を $x_j = 2\pi j/N$（$j = 0, \ldots, N-1$）とする。連続な周期関数 $u$ について

$$\tilde u_k = \frac1N\sum_{j=0}^{N-1}u(x_j)e^{-ikx_j}\quad(k \in \mathbb{Z}),\qquad I_Nu = \sum_{|k|\le K}\tilde u_k\,e_k$$

を離散 Fourier 係数、三角補間という。整数 $n$ に対し、$n \equiv n^{*} \pmod N$ かつ $|n^{*}| \le K$ を満たすただ 1 つの整数を $n^{*}$ と書き、$n \ne n^{*}$ のとき、波数 $n$ は波数 $n^{*}$ に折り返される（エイリアシング）という。

補題 3.2（指標の離散的な直交性と補間の一意性）✓

$N = 2K+1$、$x_j = 2\pi j/N$ とする。

1. 整数 $n$ について、$\frac1N\sum_{j=0}^{N-1}e^{inx_j}$ は $n$ が $N$ の倍数のとき 1、そうでないとき 0 である。
2. $v, w \in S_N$ について $\frac1N\sum_{j=0}^{N-1}v(x_j)\overline{w(x_j)} = \langle v, w\rangle$。
3. 連続な周期関数 $u$ について、$I_Nu$ は $p(x_j) = u(x_j)$（$j = 0, \ldots, N-1$）を満たすただ 1 つの $p \in S_N$ である。

証明（補題 3.2）

1. $\omega = e^{2\pi in/N}$ とおくと和は $\frac1N\sum_{j=0}^{N-1}\omega^j$ です。$n$ が $N$ の倍数なら $\omega = 1$ で和は 1。そうでなければ $\omega \ne 1$、$\omega^N = 1$ なので、等比数列の和の公式から $\sum_{j=0}^{N-1}\omega^j = (\omega^N - 1)/(\omega - 1) = 0$ です。

2. $v = \sum_{|k|\le K}a_ke_k$、$w = \sum_{|l|\le K}b_le_l$ と書くと

$$\frac1N\sum_jv(x_j)\overline{w(x_j)} = \sum_{k,l}a_k\bar b_l\cdot\frac1N\sum_je^{i(k-l)x_j}$$

です。$|k - l| \le 2K \lt N$ なので、$k - l$ が $N$ の倍数になるのは $k = l$ のときだけです。1 から右辺は $\sum_ka_k\bar b_k$ で、これは $\{e_k\}$ の正規直交性から $\langle v, w\rangle$ に等しくなります。

3. $p = I_Nu$ が補間の条件を満たすことを確かめます。

$$I_Nu(x_l) = \sum_{|k|\le K}\frac1N\sum_ju(x_j)e^{ik(x_l - x_j)} = \sum_ju(x_j)\cdot\frac1N\sum_{k=-K}^{K}\omega_{lj}^k,\qquad \omega_{lj} = e^{2\pi i(l-j)/N}$$

です。$l = j$ なら $\omega_{lj} = 1$ で内側の和は 1 です。$l \ne j$ なら $\omega_{lj} \ne 1$、$\omega_{lj}^N = 1$ で、$\sum_{k=-K}^{K}\omega_{lj}^k = \omega_{lj}^{-K}(\omega_{lj}^N - 1)/(\omega_{lj} - 1) = 0$ です。したがって $I_Nu(x_l) = u(x_l)$。一意性: $p, q \in S_N$ が同じ補間の条件を満たせば、$p - q \in S_N$ はすべての格子点で 0 なので、2 から $\|p - q\|^2 = \frac1N\sum_j|(p-q)(x_j)|^2 = 0$ です。

∎

2 は、$S_N$ の元どうしの内積が台形則で厳密に計算できることを言っています。周期関数の台形則が三角多項式に対して厳密であること（定理 3.2［数値積分 II］）の言い換えです。

次の補題は、格子点の値だけから作った係数 $\tilde u_k$ が、本来の係数 $\hat u(k)$ とどう違うかを与えます。たとえば $N = 5$ なら、波数 $3$ の波は格子点の上で波数 $-2$ の波と同じ値をとるので、$3^{*} = -2$ です。具体的な数値は §4 の例で確かめます。

補題 3.3（折り返しの公式）✓

$u$ を連続な周期関数で $\sum_n|\hat u(n)| \lt \infty$ を満たすものとする。このとき、すべての整数 $k$ について

$$\tilde u_k = \sum_{m\in\mathbb{Z}}\hat u(k + mN)$$

が成り立つ。また $I_Ne_n = e_{n^{*}}$（$n \in \mathbb{Z}$）で、

$$I_Nu = \sum_{n\in\mathbb{Z}}\hat u(n)\,e_{n^{*}},\qquad u - I_Nu = \sum_{|n| \gt K}\hat u(n)\,(e_n - e_{n^{*}})$$

である（どちらも一様に絶対収束する）。

証明（補題 3.3）

補題 2.2 から、すべての点で $u(x) = \sum_n\hat u(n)e^{inx}$ です。$x = x_j$ を代入して定義に入れ、有限和と絶対収束する級数を入れ替えると

$$\tilde u_k = \frac1N\sum_{j=0}^{N-1}\sum_n\hat u(n)e^{i(n-k)x_j} = \sum_n\hat u(n)\cdot\frac1N\sum_{j=0}^{N-1}e^{i(n-k)x_j}$$

です。補題 3.2 の 1 から、内側の和は $n - k$ が $N$ の倍数のとき 1、そうでないとき 0 なので、$n = k + mN$ の項だけが残ります。

$u = e_n$ では $\hat u(n) = 1$ で他の係数は 0 なので、$|k| \le K$ について、$\tilde u_k$ は $k \equiv n \pmod N$ すなわち $k = n^{*}$ のとき 1、そうでないとき 0 です。よって $I_Ne_n = e_{n^{*}}$ です。一般の $u$ では、最初の式から

$$I_Nu = \sum_{|k|\le K}\sum_m\hat u(k + mN)\,e_k = \sum_n\hat u(n)\,e_{n^{*}}$$

です（$n = k + mN$ と書くと $n^{*} = k$ で、$(k, m)$ と $n$ は 1 対 1 に対応する）。$|n| \le K$ では $n^{*} = n$ なので、$u = \sum_n\hat u(n)e_n$ との差は $|n| \gt K$ の項だけが残ります。$|e_n - e_{n^{*}}| \le 2$ と $\sum|\hat u(n)| \lt \infty$ から収束は一様かつ絶対です。

∎

標本化定理（定理 2.1［標本化定理とウェーブレット］）は、帯域が制限された関数は十分細かい標本から復元できると述べ、Nyquist 率より粗い標本では復元できない例がありました（命題 2.2［標本化定理とウェーブレット］）。折り返しの公式は、その周期関数版を定量的にしたものです。$N$ 点の格子が区別できる波数は $|k| \le K$ の $N$ 個だけで、それより高い波数の成分は消えるのではなく、低い波数の成分に加わります。

次の定理の証明は、誤差を「低い波数の係数に高い波数が折り返された分」と「$K$ より高い波数を捨てた分」の 2 つに分け、それぞれを Cauchy–Schwarz の不等式で抑える流れです。

定理 3.4（三角補間の誤差）✓

$s \gt 1/2$、$u \in H^s_{\mathrm{per}}$、$N = 2K+1$ とする。このとき

$$\|u - I_Nu\| \le (1 + 2Z_s)^{1/2}\Bigl(\frac2N\Bigr)^{s}|u|_{H^s}$$

が成り立つ。

証明（定理 3.4）

補題 2.2 から $\sum|\hat u(n)| \lt \infty$ なので 補題 3.3 が使えます。$I_Nu$ の Fourier 係数は、$|k| \le K$ で $\tilde u_k = \hat u(k) + \sum_{m\ne0}\hat u(k+mN)$、$|k| \gt K$ で 0 です。Parseval の等式（定理 5.1［フーリエ級数］）から

$$\|u - I_Nu\|^2 = \underbrace{\sum_{|k|\le K}\Bigl|\sum_{m\ne0}\hat u(k+mN)\Bigr|^2}_{\text{折り返しの誤差}} + \underbrace{\sum_{|n| \gt K}|\hat u(n)|^2}_{\text{切断の誤差}}$$

です。第 2 項は 定理 2.3（$r = 0$）の証明から $(2/N)^{2s}\sum_{|n| \gt K}|n|^{2s}|\hat u(n)|^2$ 以下です。

第 1 項を評価します。$|k| \le K$、$m \ne 0$ のとき $|k + mN| \ge |m|N - K \ge |m|N - N/2 \ge |m|N/2$ なので

$$\sum_{m\ne0}|k+mN|^{-2s} \le 2\sum_{m=1}^{\infty}\Bigl(\frac{mN}{2}\Bigr)^{-2s} = 2Z_s\Bigl(\frac2N\Bigr)^{2s}$$

です。Cauchy–Schwarz の不等式（命題 2.2［Hilbert 空間の基礎］）から

$$\Bigl|\sum_{m\ne0}\hat u(k+mN)\Bigr|^2 \le \sum_{m\ne0}|k+mN|^{-2s}\cdot\sum_{m\ne0}|k+mN|^{2s}|\hat u(k+mN)|^2 \le 2Z_s\Bigl(\frac2N\Bigr)^{2s}\sum_{m\ne0}|k+mN|^{2s}|\hat u(k+mN)|^2$$

です。$|k| \le K$ について足します。$|n| \gt K$ である整数 $n$ は $n = k + mN$（$|k| \le K$、$m \ne 0$）とただ 1 通りに書けるので、右辺の和は $\sum_{|n| \gt K}|n|^{2s}|\hat u(n)|^2$ になります。2 つの項を合わせて

$$\|u - I_Nu\|^2 \le (2Z_s + 1)\Bigl(\frac2N\Bigr)^{2s}\sum_{|n| \gt K}|n|^{2s}|\hat u(n)|^2 \le (1 + 2Z_s)\Bigl(\frac2N\Bigr)^{2s}|u|_{H^s}^2$$

です。

∎

補間の誤差は切断の誤差と同じ次数 $N^{-s}$ です。ただし仮定が $s \gt 1/2$ に強まっています。補間は点での値を使うので、$u$ が連続でなければ意味をもたないからです。

4. 解析関数に対する指数収束

$u$ が実軸のまわりで正則なら、Fourier 係数は指数的に減衰します（補題 3.1［数値積分 II］）。これを前節の式に入れると、誤差も指数的に減ります。

補題 4.1（等比級数の裾）✓

$0 \lt q \lt 1$、$K \ge 0$ とする。このとき

$$\sum_{n \gt K}q^n = \frac{q^{K+1}}{1-q},\qquad \sum_{n \gt K}n\,q^n \le \frac{(K+1)\,q^{K+1}}{(1-q)^2}$$

である。

証明（補題 4.1）

1 つ目は等比級数の和（命題 3.6［級数と収束判定］）に $q^{K+1}$ を掛けたものです。2 つ目を示します。まず $\sum_{j\ge1}jq^j = \sum_{i\ge1}\sum_{j\ge i}q^j = \sum_{i\ge1}q^i/(1-q) = q/(1-q)^2$ です（非負の項の二重級数なので和の順序を入れ替えてよい）。$n = K + j$（$j \ge 1$）と書くと $K + j \le (K+1)j$ なので

$$\sum_{n \gt K}n\,q^n = q^K\sum_{j\ge1}(K + j)q^j \le (K+1)\,q^K\sum_{j\ge1}jq^j = \frac{(K+1)\,q^{K+1}}{(1-q)^2}$$

です。

∎

定理 4.2（解析関数に対する指数収束）✓

$a \gt 0$ とし、$u$ は帯状領域 $\{z : |\mathrm{Im}\,z| \lt a\}$ で正則で、$u(z + 2\pi) = u(z)$、そこで $|u| \le M$ とする。$0 \lt b \lt a$、$N = 2K+1$ とすると

$$\begin{aligned} \|u - P_Nu\|_\infty &\le \frac{2M}{e^{b}-1}\,e^{-bK},\\ \|u - I_Nu\|_\infty &\le \frac{4M}{e^{b}-1}\,e^{-bK},\\ \|u' - (I_Nu)'\|_\infty &\le \frac{4M\,e^{b}}{(e^{b}-1)^2}\,(K+1)\,e^{-bK} \end{aligned}$$

が成り立つ。$K = (N-1)/2$ なので、右辺はそれぞれ $O(e^{-bN/2})$、$O(e^{-bN/2})$、$O(Ne^{-bN/2})$ である。

証明（定理 4.2）

補題 3.1［数値積分 II］ から $|\hat u(n)| \le Me^{-b|n|}$ です。$q = e^{-b}$ とおきます。とくに $\sum_n|n||\hat u(n)| \lt \infty$ です。

切断: 補題 2.2 から $u = \sum_n\hat u(n)e_n$（一様収束）なので、補題 4.1 から

$$|u(x) - P_Nu(x)| \le \sum_{|n| \gt K}|\hat u(n)| \le 2M\sum_{n \gt K}q^n = \frac{2Mq^{K+1}}{1-q} = \frac{2M}{e^b - 1}e^{-bK}$$

です。

補間: 補題 3.3 の $u - I_Nu = \sum_{|n| \gt K}\hat u(n)(e_n - e_{n^{*}})$ と $|e_n - e_{n^{*}}| \le 2$ から、切断の評価の 2 倍になります。

導関数: 級数 $\sum_n\hat u(n)e_n$ を項別に微分した級数 $\sum_nin\hat u(n)e_n$ は $\sum|n||\hat u(n)| \lt \infty$ から一様収束するので、$u$ は $C^1$ 級で $u' = \sum_nin\hat u(n)e_n$ です（定理 7.3［関数列と一様収束］ を実部と虚部に使う）。同じ理由で $\sum_n\hat u(n)e_{n^{*}}$ も項別に微分でき（$|e_{n^{*}}'| = |n^{*}| \le |n|$）、

$$u' - (I_Nu)' = \sum_{|n| \gt K}\hat u(n)\,i\,(n\,e_n - n^{*}e_{n^{*}})$$

です。$|n^{*}| \le K \lt |n|$ なので各項の絶対値は $2|n||\hat u(n)|$ 以下で、補題 4.1 から

$$\|u' - (I_Nu)'\|_\infty \le 2\sum_{|n| \gt K}|n||\hat u(n)| \le 4M\sum_{n \gt K}nq^n \le \frac{4M(K+1)q^{K+1}}{(1-q)^2} = \frac{4Me^{b}}{(e^b-1)^2}(K+1)e^{-bK}$$

です（最後は分母と分子に $e^{2b}$ を掛けた）。

∎

格子点での $(I_Nu)'$ の値が、次節のスペクトル微分行列による計算結果です。

例 4.3（スペクトル微分と中心差分）

$u(x) = e^{\sin x}$ は整関数で、$|e^{\sin z}| \le e^{\cosh a}$（$|\mathrm{Im}\,z| \le a$）なので、任意の $a$ で 定理 4.2 が使えます。一方 $u(x) = |\sin x|^3$ は $C^2$ 級ですが 3 階導関数が $x = 0, \pi$ で跳びます。それぞれの導関数を、2 次の中心差分、4 次の中心差分、スペクトル微分で近似し、格子点での最大誤差を比べます。

# Python 3.9.6 / NumPy 2.0.2
import numpy as np

def spectral_derivative(v):
    """奇数 N 点（x_j = 2πj/N）の格子値から、三角補間の導関数の格子値を返す。"""
    N = len(v)
    k = np.fft.fftfreq(N, d=1.0 / N)          # 0, 1, ..., K, -K, ..., -1
    return np.real(np.fft.ifft(1j * k * np.fft.fft(v)))

def fd2(v, h):
    return (np.roll(v, -1) - np.roll(v, 1)) / (2 * h)

def fd4(v, h):
    return (8 * (np.roll(v, -1) - np.roll(v, 1)) - (np.roll(v, -2) - np.roll(v, 2))) / (12 * h)

cases = [("exp(sin x)", lambda x: np.exp(np.sin(x)), lambda x: np.cos(x) * np.exp(np.sin(x))),
         ("|sin x|^3", lambda x: np.abs(np.sin(x)) ** 3, lambda x: 3 * np.sin(x) * np.abs(np.sin(x)) * np.cos(x))]
for name, f, df in cases:
    print(name)
    print("    N   2nd-order FD   4th-order FD     spectral")
    for N in [9, 17, 33, 65, 129, 257]:
        x = 2 * np.pi * np.arange(N) / N
        h = 2 * np.pi / N
        v, exact = f(x), df(x)
        errs = [np.max(np.abs(d - exact)) for d in (fd2(v, h), fd4(v, h), spectral_derivative(v))]
        print(f"{N:5d}   {errs[0]:.3e}      {errs[1]:.3e}      {errs[2]:.3e}")

exp(sin x)
    N   2nd-order FD   4th-order FD     spectral
    9   2.795e-01      1.186e-01      4.946e-03
   17   8.909e-02      1.318e-02      1.891e-07
   33   2.435e-02      1.033e-03      3.195e-15
   65   6.322e-03      7.133e-05      1.121e-14
  129   1.608e-03      4.642e-06      8.682e-14
  257   4.054e-04      2.952e-07      1.064e-13
|sin x|^3
    N   2nd-order FD   4th-order FD     spectral
    9   4.696e-01      2.781e-01      7.078e-02
   17   1.496e-01      7.070e-02      1.911e-02
   33   4.238e-02      1.598e-02      5.558e-03
   65   1.104e-02      3.817e-03      1.516e-03
  129   2.814e-03      9.468e-04      3.973e-04
  257   7.099e-04      2.371e-04      1.018e-04

$e^{\sin x}$ では、差分の誤差が $N$ を 2 倍にするごとに約 $1/4$、約 $1/16$ になるのに対し、スペクトル微分は $N = 33$ で丸めの大きさに達します。それ以降は誤差が少しずつ増えています。微分は波数 $k$ の成分の丸め誤差を $k$ 倍するので、丸め誤差の寄与が $N$ に比例して増えるためです。$|\sin x|^3$ では 3 つとも $N$ を 2 倍にするごとに約 $1/4$ で、方法の次数ではなく関数の滑らかさが誤差を決めています。

例 4.4（高い波数が低い波数に見える）

$N = 9$（$K = 4$）の格子で $\sin(nx)$ を標本化し、離散 Fourier 係数が 0 でない波数を調べます。補題 3.3 によれば、$\sin(nx)$ は格子点の上で $\sin(n^{*}x)$ と一致するはずです。

# Python 3.9.6 / NumPy 2.0.2
import numpy as np

N = 9                                    # 格子点の数（K = 4）
x = 2 * np.pi * np.arange(N) / N
k = np.fft.fftfreq(N, d=1.0 / N).astype(int)
for n in [3, 7, 11, 13, 22]:
    v = np.sin(n * x)
    c = np.fft.fft(v) / N                # 離散 Fourier 係数
    support = sorted(int(kk) for kk, cc in zip(k, c) if abs(cc) > 1e-12)
    alias = (n + N // 2) % N - N // 2     # n と N を法として合同な |k| <= K の整数
    diff = np.max(np.abs(v - np.sin(alias * x)))
    print(f"n = {n:2d}: nonzero at k = {support}, alias = {alias:2d}, "
          f"max_j |sin(n x_j) - sin(alias x_j)| = {diff:.1e}")

n =  3: nonzero at k = [-3, 3], alias =  3, max_j |sin(n x_j) - sin(alias x_j)| = 0.0e+00
n =  7: nonzero at k = [-2, 2], alias = -2, max_j |sin(n x_j) - sin(alias x_j)| = 4.9e-15
n = 11: nonzero at k = [-2, 2], alias =  2, max_j |sin(n x_j) - sin(alias x_j)| = 4.9e-15
n = 13: nonzero at k = [-4, 4], alias =  4, max_j |sin(n x_j) - sin(alias x_j)| = 2.4e-15
n = 22: nonzero at k = [-4, 4], alias =  4, max_j |sin(n x_j) - sin(alias x_j)| = 8.1e-15

$7 = -2 + 9$、$11 = 2 + 9$、$13 = 4 + 9$、$22 = 4 + 2\cdot9$ のとおりです。$\sin(7x)$ は格子の上では $\sin(-2x) = -\sin(2x)$ と区別できません。

5. スペクトル微分行列

定義 5.1（スペクトル微分行列）↗

$N = 2K+1$、$x_j = 2\pi j/N$ とする。格子値のベクトル $\boldsymbol{v} = (v_0, \ldots, v_{N-1})$ に対し、$p(x_j) = v_j$ を満たす $p \in S_N$（補題 3.2 の 3）をとり、$(D\boldsymbol{v})_j = p'(x_j)$ で定まる $N$ 次の行列 $D$ をスペクトル微分行列という。

命題 5.2（スペクトル微分行列の成分と固有値）✓

定義 5.1 の $D$ について次が成り立つ。

1. $D_{jj} = 0$ で、$j \ne l$ のとき $$D_{jl} = \frac{(-1)^{j-l}}{2\sin\bigl((x_j - x_l)/2\bigr)}.$$
2. $D$ は実の歪対称行列（$D^{\mathsf T} = -D$）で、巡回行列である（$D_{jl}$ は $j - l$ を $N$ で割った余りだけで決まる）。
3. $\boldsymbol{v}_k = (e^{ikx_j})_{j=0}^{N-1}$（$|k| \le K$）は互いに直交する $D$ の固有ベクトルで、$D\boldsymbol{v}_k = ik\,\boldsymbol{v}_k$。$D$ の固有値は $ik$（$|k| \le K$）ですべてであり、$\|D\|_2 = K$。
4. 2 階微分の行列（$p''(x_j)$ を与える行列）は $D^2$ に等しく、実対称で、固有値は $-k^2$（$|k| \le K$）、スペクトル半径は $K^2 = (N-1)^2/4$ である。

証明（命題 5.2）

1. $C(x) = \frac1N\sum_{|k|\le K}e^{ikx}$ とおくと、定義 3.1 から

$$p(x) = \sum_{|k|\le K}\frac1N\sum_lv_le^{ik(x - x_l)} = \sum_{l=0}^{N-1}v_l\,C(x - x_l)$$

なので、$D_{jl} = C'(x_j - x_l)$ です。Dirichlet 核の公式（補題 3.1［フーリエ級数］）から、$x$ が $2\pi$ の整数倍でないとき $C(x) = \sin(Nx/2)/(N\sin(x/2))$ で、

$$C'(x) = \frac{\frac N2\cos\frac{Nx}2\,\sin\frac x2 - \frac12\sin\frac{Nx}2\,\cos\frac x2}{N\sin^2\frac x2}$$

です。$x = x_j - x_l = 2\pi(j-l)/N$（$j \ne l$）では $Nx/2 = \pi(j-l)$ なので $\sin(Nx/2) = 0$、$\cos(Nx/2) = (-1)^{j-l}$ で、$C'(x) = (-1)^{j-l}/(2\sin(x/2))$ となります。$j = l$ では $C'(0) = \frac1N\sum_{|k|\le K}ik = 0$ です。

2. $j$ と $l$ を入れ替えると、$(-1)^{j-l}$ は変わらず、$\sin((x_j - x_l)/2)$ は符号が変わるので $D_{lj} = -D_{jl}$。$j - l$ を $j - l + N$ に置き換えると、$(-1)^{j-l}$ は $N$ が奇数なので符号が変わり、$\sin(\pi(j-l)/N + \pi) = -\sin(\pi(j-l)/N)$ も符号が変わるので、$D_{jl}$ は変わりません。

3. $e_k \in S_N$ は自分自身の補間なので、定義から $(D\boldsymbol{v}_k)_j = e_k'(x_j) = ik\,e^{ikx_j}$ です。補題 3.2 の 2 から $\frac1N\boldsymbol{v}_l^{*}\boldsymbol{v}_k = \langle e_k, e_l\rangle$ で、これは $k = l$ のとき 1、そうでないとき 0 です。したがって $\{\boldsymbol{v}_k/\sqrt N\}_{|k|\le K}$ は $\mathbb{C}^N$ の正規直交基底で、$D$ の固有値は $ik$（$|k| \le K$）の $N$ 個ですべてです。$\boldsymbol{x} = \sum_kc_k\boldsymbol{v}_k/\sqrt N$ と展開すると $\|D\boldsymbol{x}\|_2^2 = \sum_kk^2|c_k|^2 \le K^2\|\boldsymbol{x}\|_2^2$ で、$\boldsymbol{x} = \boldsymbol{v}_K$ で等号が成り立つので $\|D\|_2 = K$ です。

4. $p \in S_N$ なら $p' \in S_N$ なので、$p'$ の格子値 $D\boldsymbol{v}$ を補間する $S_N$ の元は $p'$ 自身で、その導関数 $p''$ の格子値は $D(D\boldsymbol{v})$ です。$D^2$ は $(D^2)^{\mathsf T} = (D^{\mathsf T})^2 = D^2$ から対称で、3 から $D^2\boldsymbol{v}_k = -k^2\boldsymbol{v}_k$ です。絶対値が最大の固有値は $-K^2$ です。

∎

$D$ は密行列ですが、3 から、$D\boldsymbol{v}$ は「離散 Fourier 変換をして $ik$ を掛けて逆変換する」ことで計算できます。例 4.3 の spectral_derivative がそれです。

4 は時間発展の方程式を陽的に解くときの刻み幅を決めます。熱方程式 $u_t = u_{xx}$ を空間についてスペクトル法で離散化した系 $\boldsymbol{v}' = D^2\boldsymbol{v}$ に陽的 Euler 法を使うと、波数 $k$ の成分は 1 ステップで $1 - \Delta t\,k^2$ 倍されるので、すべての成分が増えないためには $\Delta t \le 2/K^2 = 8/(N-1)^2$ が必要十分です。2 次の中心差分の場合の条件 $\Delta t \le h^2/2 = 2\pi^2/N^2$ と同じ $O(N^{-2})$ ですが、定数は $\pi^2/4 \approx 2.5$ 倍厳しくなります。差分の最大固有値 $4/h^2$ が、厳密な $K^2 = (\pi/h)^2$ より小さいためです。

注意 5.3（N が偶数のとき）

numpy.fft を使う計算では $N$ を偶数（とくに 2 の冪）にとることが多いです。$N = 2K$ のとき、波数 $K$ と $-K$ は格子の上で同じ値 $(-1)^j$ をとるので、補間の空間を $\mathrm{span}\{e_k : |k| \lt K\} \oplus \mathrm{span}\{\cos Kx\}$ にとります。このとき補間は $p(x) = \sum_lv_lC(x - x_l)$、$C(x) = \frac1N\sin(Kx)\cot(x/2)$ と書け、同じ計算で $D_{jl} = \frac12(-1)^{j-l}\cot\bigl((x_j - x_l)/2\bigr)$（$j \ne l$）が得られます。$\cos Kx$ の導関数 $-K\sin Kx$ は格子点で 0 なので、$D$ の固有値は $ik$（$|k| \lt K$）と、波数 $K$ に対応する 0 です。1 階微分では波数 $K$ の成分に 0 を掛け、2 階微分では $-K^2$ を掛けるのが整合的な扱いで、このとき 2 階微分の行列は $D^2$ とは波数 $K$ の成分だけ異なり、スペクトル半径は $N^2/4$ です。この章の定理は奇数の $N$ で述べますが、偶数でも波数 $K$ の項の扱いが変わるだけで、誤差の次数は同じです（Trefethen の教科書の第 3 章に偶数の場合の導出があります）。

6. Galerkin 法と選点法

周期境界の境界値問題

$$-u'' + c(x)\,u = f(x)\qquad(x \in \mathbb{T})$$

を考えます。この節では関数はすべて実数値とし、$c$ は連続で $0 \lt c_0 \le c(x) \le c_1$ とします。$S_N$ も実数値の元（$\hat v(-k) = \overline{\hat v(k)}$ を満たすもの）に限ります。実数値の関数の $P_Nu$ と $I_Nu$ は実数値です。

定義 6.1（Fourier–Galerkin 法と選点法）↗

$a(u, v) = \langle u', v'\rangle + \langle cu, v\rangle$（$u, v \in H^1_{\mathrm{per}}$）とおく。$f \in L^2(\mathbb{T})$ に対し、$a(u, v) = \langle f, v\rangle$（すべての $v \in H^1_{\mathrm{per}}$）を満たす $u \in H^1_{\mathrm{per}}$ を弱解という。

• $u_N \in S_N$ で $a(u_N, v) = \langle f, v\rangle$（すべての $v \in S_N$）を満たすものを Fourier–Galerkin 解という。
• $f$ が連続のとき、$u_N^{c} \in S_N$ で $-(u_N^{c})''(x_j) + c(x_j)u_N^{c}(x_j) = f(x_j)$（$j = 0, \ldots, N-1$）を満たすものを選点法の解（擬スペクトル法の解）という。

$\alpha = \min(1, c_0)$、$M_a = \max(1, c_1)$ とおくと、$\|v\|_{H^1}^2 = \|v\|^2 + \|v'\|^2$ から

$$a(v, v) \ge \alpha\|v\|_{H^1}^2,\qquad |a(u, v)| \le M_a\|u\|_{H^1}\|v\|_{H^1}$$

です。$H^1_{\mathrm{per}}$ は重み $(1 + k^2)$ つきの $\ell^2$ 空間と等長なので Hilbert 空間で、弱解は Lax–Milgram の定理（定理 3.2［楕円型方程式の弱解と変分法］）からただ 1 つ存在します。Galerkin 解も同じ定理を $S_N$ に使ってただ 1 つ存在します（定理 2.3［有限要素法 I］）。

定理 6.2（Fourier–Galerkin 法の収束）✓

$f \in L^2(\mathbb{T})$、$u$ を弱解、$u_N$ を Fourier–Galerkin 解とする。

1. $c$ が正の定数なら $\hat u(k) = \hat f(k)/(k^2 + c)$ で、$u_N = P_Nu$ である。したがって誤差は 定理 2.3 のとおりで、$u$ が 定理 4.2 の仮定を満たせば指数的に減る。
2. 一般の $c$ では、$s \ge 1$ で $u \in H^s_{\mathrm{per}}$ のとき $$\|u - u_N\|_{H^1} \le \frac{M_a}{\alpha}\|u - P_Nu\|_{H^1} \le \frac{\sqrt2\,M_a}{\alpha}(K+1)^{1-s}|u|_{H^s}.$$

証明（定理 6.2）

1. $c$ が定数のとき、弱解の定義で $v$ として $e_k$ の実部と虚部をとると、Parseval の等式から $(k^2 + c)\hat u(k) = \hat f(k)$ です。Galerkin 解についても、$|k| \le K$ で同じ式 $(k^2 + c)\hat u_N(k) = \hat f(k)$ が成り立ち、$|k| \gt K$ では $\hat u_N(k) = 0$ です。よって $u_N = P_Nu$ です。

2. Céa の補題（定理 2.5［有限要素法 I］）から $\|u - u_N\|_{H^1} \le (M_a/\alpha)\inf_{v\in S_N}\|u - v\|_{H^1}$ です。$v = P_Nu$ をとり、定理 2.3（$r = 1$）を使います。

∎

有限要素法では、近似空間が区分多項式なので $\inf_v\|u - v\|_{H^1}$ は $O(h^k)$ で頭打ちでした。同じ Céa の補題でも、近似空間を $S_N$ にとると、右辺は $u$ の滑らかさの分だけ小さくなります。

選点法は、積分を計算せずに格子点の値だけで方程式を立てられるので、変数係数や非線形の問題で使いやすい方法です。その代わり、積 $c\,u_N$ と右辺 $f$ を格子点で標本化するときに折り返しが入ります。

次の定理は、選点法の解 $u_N^{c}$ が Galerkin 解 $u_N$ からどれだけずれるかを調べます。証明の方針は 3 段階です。まず補間の誤差が小さい（定数 $N_0$ は、$N \ge N_0$ ならその誤差が $\alpha/2$ 倍以下になるように選んだ値）ことで選点法の解の存在と一意性を示し、次に両者の差が満たす式を作り、最後に $u''$ の補間誤差でずれの大きさを評価します。

定理 6.3（選点法の解と Galerkin 解の差）✓

$c$ を $C^1$ 級、$f$ を連続とし、$\gamma = \max(\|c\|_\infty, \|c'\|_\infty)$、$N_0 = 4\sqrt{10}\,\gamma/\alpha$ とおく。$N = 2K + 1 \ge N_0$ のとき、選点法の解 $u_N^{c}$ はただ 1 つ存在し、Galerkin 解 $u_N$ との差は

$$\alpha\,\|u_N^{c} - u_N\|_{H^1} \le \bigl\|(1 - I_N)(cu_N^{c} - f)\bigr\|$$

を満たす（右辺は $cu_N^{c}$ と $f$ の折り返しを含む補間誤差である）。さらに弱解 $u$ が $s \gt 5/2$ について $H^s_{\mathrm{per}}$ に属するなら

$$\|u - u_N^{c}\|_{H^1} \le 2\|u - u_N\|_{H^1} + \frac{2}{\alpha}(1 + 2Z_{s-2})^{1/2}\Bigl(\frac2N\Bigr)^{s-2}|u|_{H^s}$$

である。

証明（定理 6.3）

選点法の方程式の書き換え。 $w \in S_N$ について、$-w'' + I_N(cw) - I_Nf$ は $S_N$ の元で、格子点での値は $-w''(x_j) + c(x_j)w(x_j) - f(x_j)$ です。補題 3.2 の 3（一意性）から、$w$ が選点法の解であることと、$S_N$ の元として $-w'' + I_N(cw) = I_Nf$ が成り立つことは同値です。

補間誤差の評価。 $w \in H^1_{\mathrm{per}}$ が $C^1$ 級なら、$(cw)' = c'w + cw'$ から $|cw|_{H^1} = \|(cw)'\| \le \gamma(\|w\| + \|w'\|) \le \sqrt2\,\gamma\|w\|_{H^1}$ です。定理 3.4（$s = 1$。$Z_1 \le 2$ なので $(1 + 2Z_1)^{1/2} \le \sqrt5$）から

$$\|(1 - I_N)(cw)\| \le \sqrt5\cdot\frac2N\cdot\sqrt2\,\gamma\,\|w\|_{H^1} = \frac{2\sqrt{10}\,\gamma}{N}\|w\|_{H^1} \le \frac\alpha2\|w\|_{H^1}\qquad(N \ge N_0)$$

です。

一意性と存在。 $f$ の格子値がすべて 0 のときの解 $w \in S_N$ は $-w'' + I_N(cw) = 0$ を満たします。$w$ との内積をとり、$\langle -w'', w\rangle = \|w'\|^2$（Parseval の等式）を使うと

$$a(w, w) = \|w'\|^2 + \langle cw, w\rangle = \langle cw - I_N(cw), w\rangle \le \frac\alpha2\|w\|_{H^1}\|w\| \le \frac\alpha2\|w\|_{H^1}^2$$

です。$a(w, w) \ge \alpha\|w\|_{H^1}^2$ と合わせて $w = 0$。選点法の方程式は未知数 $N$ 個・方程式 $N$ 個の連立 1 次方程式なので、一意性から解の存在が従います。

Galerkin 解との差。 $e = u_N^{c} - u_N \in S_N$ とおきます。$v \in S_N$ について、Galerkin 解の定義から $a(u_N, v) = \langle f, v\rangle$、選点法の方程式から $\langle -(u_N^{c})'', v\rangle = \langle I_Nf - I_N(cu_N^{c}), v\rangle$ なので

$$a(e, v) = \langle -(u_N^{c})'' + cu_N^{c}, v\rangle - \langle f, v\rangle = \bigl\langle (1 - I_N)(cu_N^{c} - f), v\bigr\rangle$$

です。$v = e$ とすると $\alpha\|e\|_{H^1}^2 \le a(e, e) \le \|(1 - I_N)(cu_N^{c} - f)\|\,\|e\|_{H^1}$ で、最初の不等式が出ます。

誤差の次数。 $u \in H^s_{\mathrm{per}}$（$s \gt 5/2$）とします。$\sum_k|k||\hat u(k)| \lt \infty$（補題 2.2 を $u'$ に使う）から $u$ は $C^1$ 級です。弱解の定義で $v$ に $e_k$ の実部と虚部をとると $k^2\hat u(k) = \widehat{(f - cu)}(k)$ なので、$u'' = cu - f$（$u''$ は Fourier 係数が $-k^2\hat u(k)$ の $L^2$ の元）で、$u'' \in H^{s-2}_{\mathrm{per}}$、$|u''|_{H^{s-2}} = |u|_{H^s}$ です。そこで

$$(1 - I_N)(cu_N^{c} - f) = (1 - I_N)u'' + (1 - I_N)\bigl(c\,(u_N^{c} - u)\bigr)$$

と分けます。第 1 項は 定理 3.4（$s - 2 \gt 1/2$）から $(1 + 2Z_{s-2})^{1/2}(2/N)^{s-2}|u|_{H^s}$ 以下、第 2 項は上の補間誤差の評価から $(\alpha/2)\|u_N^{c} - u\|_{H^1} \le (\alpha/2)(\|e\|_{H^1} + \|u_N - u\|_{H^1})$ 以下です。最初の不等式に入れて $(\alpha/2)\|e\|_{H^1}$ を左辺に移すと

$$\|e\|_{H^1} \le \frac2\alpha(1 + 2Z_{s-2})^{1/2}\Bigl(\frac2N\Bigr)^{s-2}|u|_{H^s} + \|u - u_N\|_{H^1}$$

です。$\|u - u_N^{c}\|_{H^1} \le \|u - u_N\|_{H^1} + \|e\|_{H^1}$ から結論が出ます。

∎

注意 6.4（選点法の次数について）

この評価では、選点法の $H^1$ 誤差は $O(N^{2-s})$ で、Galerkin 法の $O(N^{1-s})$（定理 6.2）より $N$ の冪が 1 つ悪くなります。$|u|_{H^s}$ だけを使う評価としては、この冪は改良できません。$c$ が定数のとき、選点法の解の Fourier 係数は $\tilde f_k/(k^2 + c) = \hat u(k) + \sum_{m\ne0}\hat f(k + mN)/(k^2 + c)$ で、高い波数の $\hat f(k + mN) = ((k+mN)^2 + c)\,\hat u(k+mN)$ が低い波数 $k$ にそのまま加わります。Galerkin 法では、同じ成分は $(k+mN)^2 + c$ で割られた形で誤差に入ります。実際、$N$ ごとに $u = 2N^{-s}\cos Nx$（$|u|_{H^s} = \sqrt2$）をとると、$f = (N^2 + c)u$ の格子値はすべて $2(N^2+c)N^{-s}$ なので選点法の解は定数 $2(N^2+c)N^{-s}/c$ で、$\|u - u_N^{c}\|_{H^1} \ge 2N^{2-s}/c$ です。一方 $u_N = P_Nu = 0$ で $\|u - u_N\|_{H^1} = \sqrt2\,(1+N^2)^{1/2}N^{-s}$ です。ただし、これは $H^s$ の単位球の上での最悪の場合の比較です。1 つの $u$ を固定したとき、Fourier 係数が $|n|^{-\sigma}$ のように単調に減る関数では、両者の差は $N^{1/2}$ 倍程度にとどまります。どちらも $u$ の滑らかさに応じて任意の次数で減り、$u$ と $f$ が解析的なら指数的に減るので、滑らかな解ではこの差は実用上問題になりません。

非線形の方程式では、折り返しは安定性にも影響します。Burgers 方程式 $u_t + (u^2/2)_x = \nu u_{xx}$ の $u^2$ のような 2 次の項は、次の方法で折り返しなしに計算できます。

定義 6.5（3/2 則による折り返しの除去）↗

$u, v \in S_N$（$N = 2K+1$）の積の切断 $P_N(uv)$ を、$M \ge 3K + 1$ 点の格子で $u$ と $v$ の値を求めて掛け、その離散 Fourier 係数のうち $|k| \le K$ のものだけを残して計算する方法を、3/2 則（零詰めによる折り返しの除去）という。$M \approx (3/2)N$ である。

命題 6.6（3/2 則は Galerkin の積を厳密に与える）✓

$u, v \in S_N$、$M \ge 3K + 1$ とし、$M$ 点の格子 $y_j = 2\pi j/M$ での離散 Fourier 係数を $\tilde w^{(M)}_k = \frac1M\sum_{j=0}^{M-1}w(y_j)e^{-iky_j}$ と書く。このとき $w = uv$ について $\tilde w^{(M)}_k = \hat w(k)$（$|k| \le K$）である。$M = N$ では一般に $\tilde w^{(N)}_k \ne \hat w(k)$ である。

証明（命題 6.6）

$w = uv$ は次数 $2K$ 以下の三角多項式なので、$|n| \gt 2K$ で $\hat w(n) = 0$ です。折り返しの公式（補題 3.3 の最初の式。証明は格子点の数が偶数でもそのまま通る）から $\tilde w^{(M)}_k = \sum_m\hat w(k + mM)$ です。$|k| \le K$、$m \ne 0$ のとき $|k + mM| \ge M - K \ge 2K + 1$ なので、$m \ne 0$ の項はすべて 0 です。$M = N$ の反例は $u = v = e_K + e_{-K}$ で、$w = e_{2K} + 2 + e_{-2K}$、$\hat w(-1) = 0$ ですが、$2K = -1 + N$ なので $\tilde w^{(N)}_{-1} = \hat w(2K) = 1$ です。

∎

例 6.7（積の係数の折り返し）

$K = 8$ の実三角多項式 $u$（係数は乱数）について、$u^2$ の Fourier 係数（$|k| \le K$）を $M$ 点の格子で求め、係数の畳み込みで求めた厳密な値と比べます。

# Python 3.9.6 / NumPy 2.0.2
import numpy as np

K = 8
rng = np.random.default_rng(0)
c = rng.normal(size=K + 1) + 1j * rng.normal(size=K + 1)
c[0] = c[0].real

def coeffs(M):
    """次数 K の実三角多項式 u の u^2 の Fourier 係数（|k| <= K）を、M 点の格子上の積から求める。"""
    x = 2 * np.pi * np.arange(M) / M
    u = c[0].real + 2 * np.real(sum(c[n] * np.exp(1j * n * x) for n in range(1, K + 1)))
    w = np.fft.fft(u * u) / M
    return np.array([w[k % M] for k in range(-K, K + 1)])

full = np.zeros(2 * K + 1, dtype=complex)            # u の係数（k = -K, ..., K）
full[K:] = c
full[:K] = np.conj(c[:0:-1])
exact = np.convolve(full, full)[K:3 * K + 1]         # 畳み込みで求めた u^2 の係数（|k| <= K）
for M in [2 * K + 1, 3 * K, 3 * K + 1, 4 * K]:
    print(f"M = {M:2d}: max error of the coefficients of u^2 = {np.max(np.abs(coeffs(M) - exact)):.2e}")

M = 17: max error of the coefficients of u^2 = 8.51e+00
M = 24: max error of the coefficients of u^2 = 5.95e-01
M = 25: max error of the coefficients of u^2 = 3.58e-14
M = 32: max error of the coefficients of u^2 = 2.49e-14

$M = 3K = 24$ ではまだ折り返しが残り、$M = 3K + 1 = 25$ で丸めの大きさになります。命題 6.6 の条件 $M \ge 3K + 1$ のとおりです。

例 6.8（Fourier 選点法と 4 次 Runge–Kutta 法）

移流方程式 $u_t + u_x = 0$ と熱方程式 $u_t = u_{xx}$ を、初期値 $u_0 = e^{\sin x}$、空間は Fourier 選点法、時間は古典的な 4 次 Runge–Kutta 法で解きます。離散 Fourier 変換をした側では、波数 $k$ の成分は $w' = \lambda_kw$（移流は $\lambda_k = -ik$、熱は $\lambda_k = -k^2$）で、Runge–Kutta 法の 1 ステップは $R(\Delta t\lambda_k)$ 倍、$R(z) = 1 + z + z^2/2 + z^3/6 + z^4/24$ です。純虚数 $z = iy$ では $|R(iy)|^2 = 1 - y^6/72 + y^8/576$ なので $|R(iy)| \le 1 \iff |y| \le 2\sqrt2$、負の実数では $R(-2.7) \approx 0.879$、$R(-3.2) \approx 1.828$ です。

# Python 3.9.6 / NumPy 2.0.2
import numpy as np

def rk4(L, v, dt, steps):
    """v' = L v（L は Fourier 側の対角成分）を古典的な 4 次 Runge–Kutta 法で進める。"""
    w = np.fft.fft(v)
    for _ in range(steps):
        k1 = L * w
        k2 = L * (w + 0.5 * dt * k1)
        k3 = L * (w + 0.5 * dt * k2)
        k4 = L * (w + dt * k3)
        w = w + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6
    return np.real(np.fft.ifft(w))

u0 = lambda x: np.exp(np.sin(x))
T = 1.0
print("advection u_t + u_x = 0, error at T = 1")
print("    N \\ steps      50          100          200          400")
for N in [9, 17, 33, 65]:
    x = 2 * np.pi * np.arange(N) / N
    k = np.fft.fftfreq(N, d=1.0 / N)
    row = [np.max(np.abs(rk4(-1j * k, u0(x), T / s, s) - u0(x - T))) for s in (50, 100, 200, 400)]
    print(f"{N:5d}        " + "   ".join(f"{e:.3e}" for e in row))

print("heat u_t = u_xx, error at T = 1 against the exact solution of the semi-discrete system")
print("    N    steps (dt K^2 = 2.7)   error        steps (dt K^2 = 3.2)   error")
for N in [17, 33, 65]:
    x = 2 * np.pi * np.arange(N) / N
    k = np.fft.fftfreq(N, d=1.0 / N)
    K2 = np.max(k**2)
    ref = np.real(np.fft.ifft(np.exp(-k**2 * T) * np.fft.fft(u0(x))))
    out = []
    for z in (2.7, 3.2):
        steps = int(np.ceil(T * K2 / z))
        out.append((steps, np.max(np.abs(rk4(-k**2, u0(x), T / steps, steps) - ref))))
    print(f"{N:5d}    {out[0][0]:6d}              {out[0][1]:.3e}      {out[1][0]:6d}              {out[1][1]:.3e}")

advection u_t + u_x = 0, error at T = 1
    N \ steps      50          100          200          400
    9        1.086e-03   1.086e-03   1.086e-03   1.086e-03
   17        4.524e-08   1.836e-08   1.774e-08   1.771e-08
   33        3.304e-08   2.067e-09   1.292e-10   8.078e-12
   65        3.291e-08   2.052e-09   1.284e-10   8.033e-12
heat u_t = u_xx, error at T = 1 against the exact solution of the semi-discrete system
    N    steps (dt K^2 = 2.7)   error        steps (dt K^2 = 3.2)   error
   17        24              1.680e-07          20              3.448e-02
   33        95              6.155e-10          80              4.289e+04
   65       380              2.349e-12         320              2.629e+67

移流方程式では、$N = 9$ の行はステップ数を変えても誤差が変わりません（空間の誤差が支配）。$N = 33, 65$ の行は $N$ によらず、ステップ数を 2 倍にするごとに誤差が約 $1/16$ になっています（時間の誤差が支配）。$N = 33$ で空間の誤差は $10^{-12}$ より小さいので、精度を決めているのは時間積分です。熱方程式では、$\Delta tK^2 = 2.7$ では安定で、$3.2$ では最高波数の成分が毎ステップ約 $1.83$ 倍に増えて発散します。$K^2 \approx N^2/4$ なので、安定な刻み幅は $\Delta t \approx 11/N^2$ 以下です。

7. Chebyshev 選点法

周期的でない関数を三角多項式で近似すると、端での不一致が不連続として働き、収束が遅くなります（§8）。区間 $[-1,1]$ 上の周期的でない問題では、代数多項式による補間を使います。等間隔の点での高次の補間は不安定なので、端に密な点をとります。

定義 7.1（Chebyshev–Gauss–Lobatto 点と Chebyshev 微分行列）↗

$N \ge 1$ とする（この節の $N$ は多項式の次数で、偶奇は問わない）。$x_j = \cos(j\pi/N)$（$j = 0, \ldots, N$）を Chebyshev–Gauss–Lobatto 点という。$[-1,1]$ 上の連続関数 $u$ について、$p_N(x_j) = u(x_j)$（$j = 0, \ldots, N$）を満たす $N$ 次以下の多項式 $p_N$（定理 2.2［多項式補間］）を $u$ の Chebyshev 補間という。格子値のベクトル $\boldsymbol{v} = (v_0, \ldots, v_N)$ を補間する $N$ 次以下の多項式 $p$ について $(D_N\boldsymbol{v})_j = p'(x_j)$ で定まる $(N+1)$ 次の行列 $D_N$ を Chebyshev 微分行列という。境界値問題を、$p_N$ が内部の点 $x_1, \ldots, x_{N-1}$ で微分方程式を満たし、$x_0, x_N$ で境界条件を満たすように解く方法を Chebyshev 選点法という。

$x_0 = 1$、$x_N = -1$ で、点は右から左へ並びます。Chebyshev 多項式 $T_k$（定義 2.1［誤差保証つきの求積と関数近似］）は $T_k(\cos\theta) = \cos k\theta$ を満たすので（補題 6.1［多項式補間］）、$x = \cos\theta$ と置くと $T_k(x_j) = \cos(kj\pi/N)$ です。$x_j$ は、$\theta$ について等間隔な点 $\theta_j = j\pi/N$ の像で、Chebyshev 補間は偶関数 $u(\cos\theta)$ の余弦による三角補間です。そのため Fourier の場合と同じ折り返しが起こります。

補題 7.2（Chebyshev 多項式の折り返し）✓

整数 $k \ge 0$ に対し、$k$ を $2N$ で割った余りを $r$ とし、$k^{\sharp} = r$（$r \le N$ のとき）、$k^{\sharp} = 2N - r$（$r \gt N$ のとき）とおく。このとき $T_k(x_j) = T_{k^{\sharp}}(x_j)$（$j = 0, \ldots, N$）で、$T_k$ の Chebyshev 補間は $T_{k^{\sharp}}$ である。

証明（補題 7.2）

$k = 2Nm + r$ と書くと、$T_k(x_j) = \cos(kj\pi/N) = \cos(rj\pi/N + 2\pi mj) = \cos(rj\pi/N)$ です。$r \le N$ ならこれは $T_r(x_j)$。$r \gt N$ なら $\cos((2N - r)j\pi/N) = \cos(2\pi j - rj\pi/N) = \cos(rj\pi/N)$ なので $T_{2N-r}(x_j)$ に等しくなります。$0 \le k^{\sharp} \le N$ なので $T_{k^{\sharp}}$ は $N$ 次以下で、補間多項式の一意性（定理 2.2［多項式補間］）から $T_k$ の補間です。

∎

補題 7.3（Chebyshev 多項式の導関数）✓

$k \ge 0$ について、$[-1,1]$ 上で $|T_k'(x)| \le k^2$ であり、$T_k'(1) = k^2$、$T_k''(1) = k^2(k^2-1)/3$ である。また $0 \lt \theta \lt \pi$ で $T_k'(\cos\theta) = k\sin k\theta/\sin\theta$ である。

証明（補題 7.3）

$T_k(\cos\theta) = \cos k\theta$ を $\theta$ で微分すると $-\sin\theta\,T_k'(\cos\theta) = -k\sin k\theta$ なので、$0 \lt \theta \lt \pi$ で $T_k'(\cos\theta) = k\sin k\theta/\sin\theta$ です。$|\sin k\theta| \le k|\sin\theta|$ を $k$ についての帰納法で示します。$k = 0, 1$ では等号。$|\sin(k+1)\theta| = |\sin k\theta\cos\theta + \cos k\theta\sin\theta| \le |\sin k\theta| + |\sin\theta| \le (k+1)|\sin\theta|$ です。よって $(-1,1)$ で $|T_k'| \le k^2$ で、$T_k'$ は連続なので $[-1,1]$ でも成り立ちます。$\theta \to 0$ とすると $\sin k\theta/\sin\theta \to k$ なので $T_k'(1) = k^2$ です。

$g(\theta) = T_k(\cos\theta) = \cos k\theta$ は $g'' = -k^2g$ を満たします。$g' = -\sin\theta\,T_k'$、$g'' = \sin^2\theta\,T_k'' - \cos\theta\,T_k'$（$T_k$ の引数は $\cos\theta$）なので、$x = \cos\theta$ について $(1 - x^2)T_k''(x) - xT_k'(x) + k^2T_k(x) = 0$ が $(-1,1)$ で成り立ち、両辺は多項式なのですべての $x$ で成り立ちます。これを微分すると $(1-x^2)T_k''' - 3xT_k'' + (k^2 - 1)T_k' = 0$ で、$x = 1$ とすると $T_k''(1) = (k^2-1)T_k'(1)/3 = k^2(k^2-1)/3$ です。

∎

定理 7.4（Chebyshev 補間の誤差）✓

$\rho \gt 1$ とし、$u$ は Bernstein 楕円 $E_\rho$ で正則で $[-1,1]$ 上で実数値とする。$1 \lt r \lt \rho$、$M_r = \max_{\partial E_r}|u|$、$C_r = \sum_{j\ge1}(1+j)^2r^{-j}$ とおく。$u$ の Chebyshev 補間 $p_N$ について

$$\max_{[-1,1]}|u - p_N| \le \frac{4M_r}{r - 1}\,r^{-N},\qquad \max_{[-1,1]}|u' - p_N'| \le 4M_rC_r\,N^2r^{-N}$$

が成り立つ。

証明（定理 7.4）

定理 3.1［誤差保証つきの求積と関数近似］ から、$u = \sum_{k\ge0}a_kT_k$（$[-1,1]$ で絶対かつ一様に収束）、$|a_k| \le 2M_rr^{-k}$（$k \ge 1$）です。

補間の表示。 $q = \sum_{k\ge0}a_kT_{k^{\sharp}}$ とおきます。$\sum_k|a_k| \lt \infty$ なので、$l = 0, \ldots, N$ ごとに $k^{\sharp} = l$ となる $k$ についての和 $\sum a_k$ は絶対収束し、$q$ は $T_0, \ldots, T_N$ の 1 次結合、すなわち $N$ 次以下の多項式です。補題 7.2 から $q(x_j) = \sum_ka_kT_k(x_j) = u(x_j)$ なので、一意性から $q = p_N$ です。$k \le N$ では $k^{\sharp} = k$ なので

$$u - p_N = \sum_{k \gt N}a_k\,(T_k - T_{k^{\sharp}})$$

です。

関数値。 $[-1,1]$ で $|T_k| \le 1$（補題 6.1［多項式補間］）なので、補題 4.1（$q = 1/r$）から

$$|u - p_N| \le 2\sum_{k \gt N}|a_k| \le 4M_r\sum_{k \gt N}r^{-k} = \frac{4M_r\,r^{-N-1}}{1 - r^{-1}} = \frac{4M_r}{r-1}r^{-N}$$

です。

導関数。 補題 7.3 から $|a_kT_k'| \le 2M_rk^2r^{-k}$ で、$\sum_kk^2r^{-k}$ は比判定法（定理 6.2［級数と収束判定］）で収束するので、$\sum_ka_kT_k'$ は $[-1,1]$ で一様収束します。項別微分（定理 7.3［関数列と一様収束］）から $u' = \sum_ka_kT_k'$ で、$p_N' = \sum_ka_kT_{k^{\sharp}}'$ も同様です（$|T_{k^{\sharp}}'| \le (k^{\sharp})^2 \le k^2$）。$k \gt N$ では $(k^{\sharp})^2 \le N^2 \lt k^2$ なので

$$|u' - p_N'| \le \sum_{k \gt N}|a_k|\bigl(k^2 + (k^{\sharp})^2\bigr) \le 4M_r\sum_{k \gt N}k^2r^{-k}$$

です。$k = N + j$（$j \ge 1$）と書くと $N + j \le N(1 + j)$ なので $\sum_{k \gt N}k^2r^{-k} \le N^2r^{-N}\sum_{j\ge1}(1+j)^2r^{-j} = C_rN^2r^{-N}$ です。

∎

Chebyshev 多項式の零点（Chebyshev 点）での補間については、Lebesgue 定数を使った評価 $O(r^{-n}\log n)$ がありました（系 6.7［多項式補間］）。ここでは折り返しを使って、Lebesgue 定数を経由せずに対数の因子のない評価を得ています。導関数の誤差は $N^2$ 倍になります。この $N^2$ は $T_k'$ が端で $k^2$ の大きさをもつことから来ていて、次の命題のとおり微分行列の大きさにも現れます。

次の命題の 1 の証明では、点 $x_0, \ldots, x_N$ をちょうど零点にもつ多項式 $\omega(x) = (1 - x^2)T_N'(x)$ を使います。Lagrange の基底多項式は「その点で 1、他の点で 0」の多項式で、分母にこの $\omega$ の導関数が現れるため、$\omega'$ の値さえ分かれば微分行列の成分が書き下せるからです。

命題 7.5（Chebyshev 微分行列）✓

$N \ge 1$、$c_0 = c_N = 2$、$c_j = 1$（$0 \lt j \lt N$）とする。

1. $j \ne l$ のとき $D_N$ の成分は $$(D_N)_{jl} = \frac{c_j}{c_l}\cdot\frac{(-1)^{j+l}}{x_j - x_l}$$ で、対角成分は $(D_N)_{jj} = -\sum_{l\ne j}(D_N)_{jl}$ である。
2. $N$ 次以下の多項式 $p$ は $p(\cos\theta) = \sum_{k=0}^{N}b_k\cos k\theta$ と書け、$0 \lt \theta \lt \pi$ で $p'(\cos\theta) = \sum_{k=1}^{N}k\,b_k\sin k\theta/\sin\theta$ である。
3. 端の点の間隔は $x_0 - x_1 = 2\sin^2(\pi/(2N)) \le \pi^2/(2N^2)$ である。
4. $\|D_N\|_\infty \ge N^2$、$\|D_N^2\|_\infty \ge N^2(N^2-1)/3$ である。一方 $D_N^{N+1} = 0$ である。

証明（命題 7.5）

1. $\omega(x) = (1 - x^2)T_N'(x)$ とおきます。補題 7.3 から $T_N'(\cos\theta) = N\sin N\theta/\sin\theta$ で、これは $\theta_j = j\pi/N$（$j = 1, \ldots, N-1$）で 0 です。$T_N'$ は $N - 1$ 次なので、零点は $x_1, \ldots, x_{N-1}$ ですべてで、$\omega$ は $x_0, \ldots, x_N$ をちょうど零点とする $N + 1$ 次の多項式です。$\omega(x) = \kappa\prod_m(x - x_m)$ と書くと $\omega'(x_l) = \kappa\prod_{m\ne l}(x_l - x_m)$ なので、Lagrange の基底多項式（定理 2.2［多項式補間］）は $\ell_l(x) = \omega(x)/(\omega'(x_l)(x - x_l))$ です。$p = \sum_lv_l\ell_l$ から $(D_N)_{jl} = \ell_l'(x_j)$ で、$j \ne l$ のとき $\omega(x_j) = 0$ を使うと

$$\ell_l'(x_j) = \frac{\omega'(x_j)(x_j - x_l) - \omega(x_j)}{\omega'(x_l)(x_j - x_l)^2} = \frac{\omega'(x_j)}{\omega'(x_l)(x_j - x_l)}$$

です。$\omega'(x_j)$ を求めます。$0 \lt \theta \lt \pi$ で $\omega(\cos\theta) = N\sin\theta\sin N\theta$ なので、$\theta$ で微分して

$$\omega'(\cos\theta) = -N\Bigl(\cos\theta\,\frac{\sin N\theta}{\sin\theta} + N\cos N\theta\Bigr)$$

です。$\theta = \theta_j$（$0 \lt j \lt N$）では $\sin N\theta_j = 0$、$\cos N\theta_j = (-1)^j$ なので $\omega'(x_j) = -N^2(-1)^j$。$\theta \to 0$ では $\sin N\theta/\sin\theta \to N$ なので $\omega'(1) = -2N^2$。$\theta = \pi - \varphi$、$\varphi \to 0$ では $\sin N\theta/\sin\theta = -(-1)^N\sin N\varphi/\sin\varphi \to -(-1)^NN$、$\cos\theta \to -1$、$\cos N\theta \to (-1)^N$ なので $\omega'(-1) = -2N^2(-1)^N$。まとめて $\omega'(x_j) = -c_jN^2(-1)^j$ で、$(D_N)_{jl} = c_j(-1)^j/(c_l(-1)^l(x_j - x_l))$ です。対角成分: 定数関数 1 は $N$ 次以下の多項式で導関数が 0 なので、$D_N$ の各行の成分の和は 0 です。

2. $T_0, \ldots, T_N$ は次数が $0, \ldots, N$ なので $N$ 次以下の多項式の基底で、$p = \sum_kb_kT_k$ と書けます。$T_k(\cos\theta) = \cos k\theta$ から最初の式が出て、$\theta$ で微分すると $-\sin\theta\,p'(\cos\theta) = -\sum_kkb_k\sin k\theta$ です。

3. $1 - \cos(\pi/N) = 2\sin^2(\pi/(2N))$ と $\sin t \le t$ から出ます。

4. $\boldsymbol{v} = (T_N(x_j))_j = ((-1)^j)_j$ は $\|\boldsymbol{v}\|_\infty = 1$ で、補間する多項式は $T_N$ 自身なので、$(D_N\boldsymbol{v})_0 = T_N'(1) = N^2$（補題 7.3）。よって $\|D_N\|_\infty \ge \|D_N\boldsymbol{v}\|_\infty \ge N^2$ です。$N$ 次以下の多項式の導関数は $N$ 次以下なので、$D_N^2\boldsymbol{v}$ は 2 階導関数の格子値で、$(D_N^2\boldsymbol{v})_0 = T_N''(1) = N^2(N^2-1)/3$ です。同じ理由で $D_N^{N+1}\boldsymbol{v}$ は $N + 1$ 階導関数の格子値で、$N$ 次以下の多項式では 0 です。

∎

2 が「Fourier の場合への帰着」です。Chebyshev 補間の微分は、余弦級数の係数 $b_k$ を離散余弦変換（FFT）で求め、$kb_k$ を係数とする正弦級数を作って $\sin\theta$ で割ることでも計算できます。

4 の $D_N^{N+1} = 0$ は、境界条件を課さない $D_N$ の固有値がすべて 0 であることを意味します。時間発展の安定性に関係するのは、境界条件を入れた行列です。たとえば $u(\pm1) = 0$ のもとでの 2 階微分は、$D_N^2$ から最初と最後の行と列を除いた $(N-1)$ 次の行列 $\tilde D^{(2)}_N$ で表されます。

注意 7.6（境界条件つきの微分行列の固有値（文献））

$\tilde D^{(2)}_N$ の固有値はすべて負の実数で、絶対値が最大のものは $N \to \infty$ で約 $0.048\,N^4$ のように増えます（Gottlieb–Lustman の定理と Weideman–Trefethen の解析。参考文献 [6]）。この証明は多項式の零点についての議論を要するので、この章では下の数値で確かめるにとどめます。4 の下界 $\|D_N^2\|_\infty \ge N^2(N^2-1)/3$ は、この増え方がノルムの段階ですでに現れていることを示します。上からの評価 $\max_{[-1,1]}|p'| \le N^2\max_{[-1,1]}|p|$（$p$ は $N$ 次以下の多項式）は Markov の不等式で、Cheney の教科書（参考文献 [5]）に証明があります。対角成分の閉じた形（$(D_N)_{00} = (2N^2+1)/6$ など）は Trefethen の教科書の第 6 章にあります。

したがって、熱方程式を Chebyshev 選点法と陽的な時間積分で解くと、刻み幅は $\Delta t = O(N^{-4})$ に制限されます。1 階の方程式（移流・波動）では $O(N^{-2})$ です。Fourier の場合（$O(N^{-2})$・$O(N^{-1})$）より厳しいのは、3 のとおり端の点の間隔が $O(N^{-2})$ だからで、最も細かい格子幅についての CFL 条件に当たります。陰的な方法か、端の密集を緩める変数変換を使うのが普通です。

例 7.7（Chebyshev 選点法による境界値問題）

$u'' = e^{4x}$、$u(\pm1) = 0$ の解は $u(x) = (e^{4x} - x\sinh4 - \cosh4)/16$ です。$u$ は整関数なので、定理 7.4 から補間の誤差はどんな $r$ についても $O(r^{-N})$ です。選点法の方程式 $\tilde D^{(2)}_N\boldsymbol{u} = (e^{4x_j})_{j=1}^{N-1}$ を解き、あわせて $\tilde D^{(2)}_N$ のスペクトル半径と $\|D_N\|_\infty$ を調べます。

# Python 3.9.6 / NumPy 2.0.2
import numpy as np

def cheb(N):
    """Chebyshev–Gauss–Lobatto 点 x_j = cos(jπ/N) と微分行列 D（(N+1) 次）。"""
    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
    dx = x[:, None] - x[None, :] + np.eye(N + 1)       # 対角は仮に 1
    D = (c[:, None] / c[None, :]) / dx
    D -= np.diag(D.sum(axis=1))                        # 対角成分は行和が 0 になるように決める
    return D, x

exact = lambda x: (np.exp(4 * x) - x * np.sinh(4.0) - np.cosh(4.0)) / 16
print("   N    max error    rho(D2_interior)   rho / N^4    ||D||_inf / N^2")
for N in [4, 8, 12, 16, 20, 24, 32, 48, 64]:
    D, x = cheb(N)
    D2 = (D @ D)[1:N, 1:N]                             # 両端の行と列を除く（u(±1) = 0）
    u = np.linalg.solve(D2, np.exp(4 * x[1:N]))
    err = np.max(np.abs(u - exact(x[1:N])))
    rho = np.max(np.abs(np.linalg.eigvals(D2)))
    print(f"{N:4d}    {err:.3e}    {rho:.6e}       {rho / N**4:.5f}      "
          f"{np.max(np.sum(np.abs(D), axis=1)) / N**2:.4f}")

   N    max error    rho(D2_interior)   rho / N^4    ||D||_inf / N^2
   4    1.039e-01    1.954400e+01       0.07634      1.0000
   8    2.162e-04    2.143723e+02       0.05234      1.0000
  12    1.067e-07    1.023326e+03       0.04935      1.0000
  16    1.935e-11    3.174792e+03       0.04844      1.0000
  20    2.887e-15    7.687117e+03       0.04804      1.0000
  24    7.994e-15    1.586978e+04       0.04783      1.0000
  32    3.508e-14    4.993925e+04       0.04763      1.0000
  48    5.551e-14    2.520431e+05       0.04748      1.0000
  64    1.814e-13    7.957304e+05       0.04743      1.0000

誤差は $N = 20$ で丸めの大きさに達します（$N = 16$ で $2\times10^{-11}$）。それ以降は、行列の条件数が $N^4$ で増えるため、丸め誤差の寄与が少しずつ大きくなります。スペクトル半径は $0.047\,N^4$ 程度で、注意 7.6 の値と合います。$\|D_N\|_\infty/N^2$ は表示の桁で 1 で、命題 7.5 の 4 の下界 $N^2$ がこの範囲では等号に近いことが分かります。

例 7.8（調和振動子の固有値）

$-u'' + y^2u = Eu$（$y \in \mathbb{R}$）の固有値は $E = 1, 3, 5, \ldots$ です。区間 $[-L, L]$（$L = 8$）に切って両端で $u = 0$ とし、$y = Lx$ の変数変換で Chebyshev 選点法を使うと、行列 $-\tilde D^{(2)}_N/L^2 + \mathrm{diag}(L^2x_j^2)$ の固有値が近似を与えます（例 7.7 の cheb を使い、np.linalg.eigvals で固有値を求めました）。

$N$	$E_1$	$E_2$	$E_3$	$E_4$
16	0.9054614145799	3.5335247519847	4.3172265629918	10.8740876566453
24	0.9987954905584	3.0148542022036	4.8858222558542	7.3416780603733
32	0.9999972825215	3.0000635816753	4.9991772019360	7.0053199034813
48	0.9999999999998	3.0000000000101	4.9999999997456	7.0000000039444
64	1.0000000000000	3.0000000000000	5.0000000000000	7.0000000000000

固有関数 $e^{-y^2/2}\times$（多項式）は整関数なので、収束は指数的です。$N = 16$ では点が足りず、近似固有値は真の値の上にも下にもずれています。選点法は Galerkin 法ではないので、近似固有値が上界になる保証（Rayleigh–Ritz の上界）はありません。区間を切ったことによる誤差は $e^{-L^2}$ 程度で、この表では見えません（切断の誤差は固有値問題の離散化で評価します）。無限区間を切らずに扱うには、Hermite 関数を基底にとる方法があります（Boyd の教科書の第 17 章）。

8. 滑らかでない関数

スペクトル法の高い精度は、解が滑らかであることに依存します。不連続があると、誤差は不連続点の近くだけでなく区間全体で大きくなります。

命題 8.1（不連続な関数の Fourier 切断）✓

$f(\theta) = \mathrm{sgn}\,\theta$（$-\pi \lt \theta \lt \pi$）とし、$S_nf = \sum_{|k|\le n}\hat f(k)e_k$ を部分和とする。

1. $\sup_\theta|f(\theta) - S_nf(\theta)|$ は $n \to \infty$ で 0 に収束しない。
2. すべての $n \ge 0$ で $\|f - S_nf\|^2 \ge \dfrac{4}{\pi^2(n+2)}$。
3. 不連続点から離れた点 $\theta = \pi/2$ で、すべての $m \ge 0$ について $\bigl|f(\pi/2) - S_{2m+1}f(\pi/2)\bigr| \ge \dfrac{2}{\pi(2m+5)}$。

証明（命題 8.1）

1. 定理 6.1［フーリエ級数］ から $S_{2m+1}f(\pi/(2m+2)) \to \frac2\pi\int_0^\pi\frac{\sin u}{u}du \approx 1.179$ で、$f(\pi/(2m+2)) = 1$ なので、十分大きい $m$ で $\sup_\theta|f - S_{2m+1}f| \gt 0.17$ です。$\hat f(k)$ は偶数の $k$ で 0 なので（下の 2 の計算）$S_{2m+2}f = S_{2m+1}f$ で、偶数の $n$ でも同じです。

2. $\hat f(k) = \frac1{2\pi}\int_{-\pi}^{\pi}\mathrm{sgn}\,\theta\,e^{-ik\theta}d\theta = -\frac{i}{\pi}\int_0^\pi\sin k\theta\,d\theta = -\frac{i(1 - \cos k\pi)}{\pi k}$ なので、$k$ が奇数のとき $|\hat f(k)|^2 = 4/(\pi^2k^2)$、偶数のとき 0 です。Parseval の等式から $\|f - S_nf\|^2 = \frac{8}{\pi^2}\sum_{k \gt n,\ k\text{ 奇数}}k^{-2}$ です。$n$ より大きい最小の奇数を $a$ とすると $a \le n + 2$ で、$t \mapsto (a + 2t)^{-2}$ は減少するので

$$\sum_{j\ge0}\frac1{(a+2j)^2} \ge \int_0^\infty\frac{dt}{(a+2t)^2} = \frac1{2a} \ge \frac1{2(n+2)}$$

です。

3. 定理 6.1［フーリエ級数］ の式から $S_{2m+1}f(\pi/2) = \frac4\pi\sum_{k=0}^{m}\frac{(-1)^k}{2k+1}$ です。$f$ は $\pi/2$ の近くで定数なので、定理 3.4［フーリエ級数］ から $S_nf(\pi/2) \to f(\pi/2) = 1$ です。したがって

$$f(\pi/2) - S_{2m+1}f(\pi/2) = \frac4\pi\sum_{k \gt m}\frac{(-1)^k}{2k+1} = \frac4\pi(-1)^{m+1}\sum_{j\ge0}\Bigl(\frac1{2m+3+4j} - \frac1{2m+5+4j}\Bigr)$$

です（収束する級数の隣り合う 2 項をまとめた）。各項は $\frac{2}{(2m+3+4j)(2m+5+4j)} \ge \frac{2}{(2m+5+4j)^2}$ で、2 と同じく積分と比べると $\sum_{j\ge0}\frac{2}{(a+4j)^2} \ge 2\int_0^\infty\frac{dt}{(a+4t)^2} = \frac1{2a}$（$a = 2m+5$）です。よって絶対値は $\frac4\pi\cdot\frac1{2(2m+5)}$ 以上です。

∎

3 が差分法との違いです。差分法では、不連続の影響は（線形で安定なスキームなら）不連続点の近くにとどまります。スペクトル法では、どの点の近似にもすべての点の値を使うので、1 か所の不連続が区間全体の誤差を $O(1/N)$ にします。周期的でない滑らかな関数（たとえば $[0, 2\pi]$ 上の $x$）を Fourier 級数で近似しようとした場合も、周期的に延長すると端で跳ぶので、同じことが起こります。これが、周期的でない問題に Chebyshev 多項式を使う理由です。

注意 8.2（この章で扱わなかったこと）

区間を要素に分けて各要素で高次の多項式を使うスペクトル要素法、Legendre 多項式による Galerkin 法とその安定性、不連続のある解に対するフィルタや Gegenbauer 再構成、非線形の方程式で折り返しが引き起こす不安定性の理論は、Canuto・Hussaini・Quarteroni・Zang の教科書（参考文献 [2]）と Gottlieb–Orszag のモノグラフ（参考文献 [4]）にあります。折り返しの除去は Orszag の 1971 年の論文（参考文献 [7]）が出発点です。

9. 演習

演習 9.1易

$N = 5$（$K = 2$）とする。$3^{*}$ を求め、$\sin 3x$ と $-\sin 2x$ が 5 個の格子点 $x_j = 2\pi j/5$ で一致することを確かめよ。また $\|e_3 - I_5e_3\|$ を求めよ。

解答

$3 \equiv -2 \pmod 5$ で $|-2| \le 2$ なので $3^{*} = -2$。$3x_j = -2x_j + 5x_j = -2x_j + 2\pi j$ なので $\sin 3x_j = \sin(-2x_j) = -\sin 2x_j$ です。補題 3.3 から $I_5e_3 = e_{-2}$ で、$e_3$ と $e_{-2}$ は直交する単位ベクトルなので $\|e_3 - e_{-2}\|^2 = 1 + 1 = 2$、$\|e_3 - I_5e_3\| = \sqrt2$ です。補間が $e_3$ を「捨てる」のではなく、別の波数に置き換えていることが、切断の誤差（$\|e_3 - P_5e_3\| = 1$）より大きい値に表れています。

演習 9.2標準

$u$ を連続な周期関数で $\sum_n|\hat u(n)| \lt \infty$ を満たすものとする。台形則 $T_N(u) = \frac{2\pi}{N}\sum_{j=0}^{N-1}u(x_j)$ と積分 $I(u) = \int_0^{2\pi}u\,dx$ について、$T_N(u) - I(u) = 2\pi\sum_{m\ne0}\hat u(mN)$ を折り返しの公式から導け。

解答

定義から $T_N(u) = 2\pi\tilde u_0$、$I(u) = 2\pi\hat u(0)$ です。補題 3.3（$k = 0$）から $\tilde u_0 = \sum_m\hat u(mN)$ なので、$T_N(u) - I(u) = 2\pi(\tilde u_0 - \hat u(0)) = 2\pi\sum_{m\ne0}\hat u(mN)$ です。これは 定理 3.2［数値積分 II］ の 2 の等式で、台形則の誤差は「波数 $mN$ の成分が波数 0 に折り返された分」です。$u$ が帯状領域で正則なら各項が $Me^{-b|m|N}$ 以下なので、誤差は $O(e^{-bN})$ になります。補間の誤差 $O(e^{-bN/2})$ より指数が 2 倍よいのは、補間では波数 $K \approx N/2$ のすぐ外の成分が効くのに対し、積分では波数 $N$ 以上の成分しか効かないからです。

演習 9.3標準

$N = 2K+1$、$K \ge 1$ とする。半離散の移流方程式 $\boldsymbol{v}' = -D\boldsymbol{v}$ に陽的 Euler 法 $\boldsymbol{v}^{n+1} = (I - \Delta tD)\boldsymbol{v}^n$ を使うと、どんな $\Delta t \gt 0$ でも $\|(I - \Delta tD)^n\|_2 \to \infty$（$n \to \infty$）となることを示せ。また、$T = n\Delta t$ を固定して $\Delta t \to 0$ とするとき $\|(I - \Delta tD)^n\|_2$ が有界にとどまるための、$\Delta t$ と $K$ の関係を求めよ。

解答

命題 5.2 の 3 から、$I - \Delta tD$ は正規直交な固有ベクトル $\boldsymbol{v}_k/\sqrt N$ と固有値 $1 - ik\Delta t$ をもつので、$\|(I - \Delta tD)^n\|_2 = \max_{|k|\le K}|1 - ik\Delta t|^n = (1 + K^2\Delta t^2)^{n/2}$ です。$K \ge 1$、$\Delta t \gt 0$ なら底が 1 より大きいので発散します。$n = T/\Delta t$ のとき $(1 + K^2\Delta t^2)^{T/(2\Delta t)} \le e^{TK^2\Delta t/2}$ で、また $K^2\Delta t \to \infty$ なら左辺は発散するので（$y = K^2\Delta t^2 \le 1$ のときは $\log(1 + y) \ge y/2$ から左辺は $e^{TK^2\Delta t/4}$ 以上、$y \gt 1$ のときは $2^{T/(2\Delta t)}$ 以上）、有界であるための条件は $K^2\Delta t$ が有界、すなわち $\Delta t = O(N^{-2})$ です。純虚数の固有値をもつ問題に陽的 Euler 法を使うと、1 階の方程式でも放物型と同じ厳しさの刻み幅が要ります。例 6.8 の 4 次 Runge–Kutta 法では、安定領域が虚軸の区間 $|y| \le 2\sqrt2$ を含むので $\Delta t \le 2\sqrt2/K$ で済みます。

演習 9.4難

$N \ge 1$ とする。$T_{N+1}$ の Chebyshev 補間 $p_N$ を求め、誤差が $T_{N+1}(x) - p_N(x) = -\frac2N(1 - x^2)T_N'(x)$ と書けることを示せ。これを補間誤差の公式（定理 4.1［多項式補間］）と比べよ。

解答

補題 7.2 で $k = N + 1$ とします。$N \ge 2$ なら $r = N + 1 \gt N$ なので $k^{\sharp} = 2N - (N+1) = N - 1$、$N = 1$ なら $r = 0$ で $k^{\sharp} = 0 = N - 1$。どちらも $p_N = T_{N-1}$ です。$x = \cos\theta$ とすると

$$T_{N+1}(x) - T_{N-1}(x) = \cos(N+1)\theta - \cos(N-1)\theta = -2\sin N\theta\,\sin\theta$$

で、補題 7.3 から $\sin N\theta\sin\theta = \sin^2\theta\cdot T_N'(\cos\theta)/N = (1 - x^2)T_N'(x)/N$ です。

補間誤差の公式は $u - p_N = \frac{u^{(N+1)}(\xi)}{(N+1)!}\prod_{j=0}^{N}(x - x_j)$ です。$u = T_{N+1}$ では $u^{(N+1)} = 2^N(N+1)!$（最高次の係数が $2^N$）なので、誤差は $2^N\prod_j(x - x_j)$ です。一方 $\omega(x) = (1 - x^2)T_N'(x)$ は $x_0, \ldots, x_N$ を零点にもち、最高次の係数は $-N\cdot2^{N-1}$ なので $\prod_j(x - x_j) = -\omega(x)/(N2^{N-1})$、よって誤差は $-2\omega(x)/N$ で、上の結果と一致します。$\max_{[-1,1]}|\prod_j(x - x_j)| = \max|2\sin N\theta\sin\theta|/2^N \le 2^{1-N}$ で、Chebyshev 点の節点多項式の最大値 $2^{-N}$（定理 6.2［多項式補間］）の 2 倍以内です。端点を補間点に含めても、節点多項式の大きさはほとんど変わりません。

参考文献

1. L. N. Trefethen, Spectral Methods in MATLAB, SIAM, 2000 — 第 3 章（周期格子と FFT）、第 4 章（滑らかさとスペクトル精度）、第 6 章（Chebyshev 微分行列）。
2. C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zang, Spectral Methods: Fundamentals in Single Domains, Springer, 2006 — Fourier・Chebyshev 近似の誤差評価と、Galerkin 法・選点法の解析。
3. J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed., Dover, 2001 — 折り返しと 3/2 則、無限区間の基底。
4. D. Gottlieb, S. A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, 1977.
5. E. W. Cheney, Introduction to Approximation Theory, 2nd ed., AMS Chelsea, 1982 — Markov の不等式。
6. J. A. C. Weideman, L. N. Trefethen, “The eigenvalues of second-order spectral differentiation matrices”, SIAM Journal on Numerical Analysis 25 (1988), 1279–1298.
7. S. A. Orszag, “On the elimination of aliasing in finite-difference schemes by filtering high-wavenumber components”, Journal of the Atmospheric Sciences 28 (1971), 1074.