# 多項式補間：誤差公式・Runge 現象・Chebyshev 点と Lebesgue 定数 > 補間多項式の存在と一意性、Newton 形と差分商、Rolle の定理による誤差公式と Hermite 補間を証明し、Lebesgue 定数で等間隔点の指数的な不安定さと Chebyshev 点の対数的な良さを評価する。 > https://rikai.mugen-giken.com/mathematics/numerical-analysis/polynomial-interpolation ## 0. この記事の要点 何点かでしか値がわからない関数を、多項式で置き換えて近似するにはどうすればよいか、という話です。点を通る多項式は簡単に作れますが、点をどう並べるかで結果が大きく変わり、等間隔に並べると点を増やすほどかえって悪くなる場合があります。この記事を読み終えると、補間多項式の作り方と誤差の見積もり方、そして点の置き方の良し悪しを数で測る方法がわかります。 - 相異なる $n+1$ 点での値を与えると、それを通る $n$ 次以下の多項式（**補間多項式**）がただ 1 つ存在します。証明は Vandermonde 行列式が $0$ でないことです。 - 補間多項式は**差分商**を係数とする Newton 形で書け、点を 1 つ足すたびに項を 1 つ足すだけで更新できます。差分商は $f^{(n)}(\xi)/n!$ に等しく、導関数の離散版です。 - **補間誤差**は $f(x) - p_n(x) = \dfrac{f^{(n+1)}(\xi)}{(n+1)!}\prod_j (x - x_j)$ です。証明は Rolle の定理の繰り返しです。値と導関数の両方を合わせる Hermite 補間も同じ形の誤差公式をもち、これが Gauss 求積の誤差評価の土台になります。 - 補間の安定性は **Lebesgue 定数** $\Lambda_n$ で測ります。$\|f - p_n\| \le (1 + \Lambda_n) E_n(f)$ で、$\Lambda_n$ はデータの誤差が何倍に増幅されるかを表す[条件数](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-condition-number)でもあります。 - 等間隔点では $\Lambda_n \ge 2^n/(4n^2)$ と指数的に増え、Runge の関数 $1/(1+25x^2)$ では補間が発散します。**Chebyshev 点**では $\Lambda_n \le (2/\pi)\log(2n+1) + 1 + 4/\pi$ と対数的にしか増えず、解析関数なら誤差は幾何級数的に減ります。 ## 1. 動機 関数 $f$ の値がいくつかの点でしかわからないとき、あるいは $f$ が複雑で何度も評価したくないとき、$f$ を多項式で置き換えて計算します。数値積分（多項式を積分する）、数値微分（多項式を微分する）、常微分方程式の多段法（過去の値を多項式で外挿する）、非線形方程式の割線法（2 点を通る 1 次式の零点）は、どれもこの置き換えの上に立っています。 自然な発想は「点を増やせば近似はよくなるだろう」というものです。ところが 1901 年に C. Runge が示したように、$[-1,1]$ の**等間隔点**で $1/(1+25x^2)$ を補間すると、点を増やすほど区間の端の近くで誤差が爆発します。滑らかな（実解析的な）関数でも、点の配置を誤ると補間は収束しないのです。一方、区間の端に点を密に置く **Chebyshev 点**を使うと、同じ関数の補間は幾何級数的に収束します。 この違いを説明するのが本章の主題です。誤差公式の中の積 $\prod_j(x - x_j)$ の大きさと、Lebesgue 定数という「補間という操作の条件数」の 2 つを調べると、等間隔点と Chebyshev 点の差が定量的にわかります。 ## 2. 補間多項式の存在と一意性 以下、$\mathcal{P}_n$ で実係数の $n$ 次以下の多項式全体を表します。 $n \ge 1$ とし、$x_0, \ldots, x_n$ を実数とする。$(i, j)$ 成分が $x_i^j$（$0 \le i, j \le n$）の行列 $V$ について $$ \det V = \prod_{0 \le i \lt j \le n} (x_j - x_i) $$ が成り立つ。とくに $x_0, \ldots, x_n$ が相異なれば $\det V \ne 0$ である。 $n$ についての帰納法で示します。$n = 1$ では $\det V = x_1 - x_0$ で正しいです。$n \ge 2$ とし、$n$ 個の点についての主張を仮定します。 $V$ の第 $j$ 列から第 $j-1$ 列の $x_0$ 倍を引く操作を、$j = n, n-1, \ldots, 1$ の順に行います（後ろから行うので、引く列はつねに元のままです）。この操作で行列式は変わりません（ の 3）。操作後の $(i, j)$ 成分は、$j \ge 1$ で $x_i^j - x_0 x_i^{j-1} = x_i^{j-1}(x_i - x_0)$ です。とくに第 $0$ 行は $(1, 0, \ldots, 0)$ になります。第 $0$ 行について余因子展開（）すると、$\det V$ は第 $i$ 行（$i = 1, \ldots, n$）が $(x_i - x_0)(1, x_i, \ldots, x_i^{n-1})$ の $n$ 次の行列式に等しくなります。各行から $x_i - x_0$ をくくり出すと（ の 2 を、5 により行に適用） $$ \det V = \prod_{i=1}^{n}(x_i - x_0) \cdot \det\bigl(x_i^{j}\bigr)_{1 \le i \le n,\ 0 \le j \le n-1} $$ です。右辺の行列式は $x_1, \ldots, x_n$ の Vandermonde 行列式なので、帰納法の仮定から $\prod_{1 \le i \lt j \le n}(x_j - x_i)$ に等しく、主張を得ます。 $x_0, \ldots, x_n$ を相異なる実数、$y_0, \ldots, y_n$ を実数とする。$p(x_i) = y_i$（$i = 0, \ldots, n$）を満たす $p \in \mathcal{P}_n$ がただ 1 つ存在し、 $$ p(x) = \sum_{i=0}^{n} y_i\,\ell_i(x), \qquad \ell_i(x) = \prod_{j \ne i}\frac{x - x_j}{x_i - x_j} $$ と書ける（**Lagrange 形**）。 $n = 0$ なら定数 $y_0$ がただ 1 つの解です。$n \ge 1$ とします。$p(x) = \sum_{j=0}^n a_j x^j$ とおくと、条件 $p(x_i) = y_i$ は係数ベクトル $a = (a_0, \ldots, a_n)^{\mathsf{T}}$ についての連立 1 次方程式 $Va = y$ です（$V$ は の行列）。$\det V \ne 0$ なので、正則性の判定（ の 1 ⇒ 5）からただ 1 つの解があります。$\ell_i$ は $n$ 次の多項式で $\ell_i(x_k) = 1$（$k = i$）、$0$（$k \ne i$）を満たすので、$\sum_i y_i\ell_i$ は条件を満たす $\mathcal{P}_n$ の元であり、一意性からこれが解です。 $q \in \mathcal{P}_n$ が相異なる $n + 1$ 個の点で $0$ になるならば、$q$ は零多項式である。 $q$ と零多項式は、どちらもその $n+1$ 点で値 $0$ をとる $\mathcal{P}_n$ の元です。 の一意性から両者は等しくなります。 以下、$f$ の値 $y_i = f(x_i)$ を補間する $\mathcal{P}_n$ の元を $p_n$（点を明示したいときは $p_{x_0, \ldots, x_n}$）と書き、$\omega(x) = \prod_{j=0}^{n}(x - x_j)$ とおきます。Lagrange 形は理論には便利ですが、点を 1 つ足すとすべての $\ell_i$ が変わります。次の節の Newton 形はこの欠点を解消します。 ## 3. 差分商と Newton 形 相異なる点 $x_0, \ldots, x_k$ と関数 $f$ について、$f$ をこの $k+1$ 点で補間する $\mathcal{P}_k$ の元の $x^k$ の係数を $f[x_0, \ldots, x_k]$ と書き、$f$ の **$k$ 階の差分商**という。$f[x_0] = f(x_0)$ である。 $f[x_0, \ldots, x_k]$ は点の並べ方によらず、 $$ f[x_0, \ldots, x_k] = \sum_{i=0}^{k}\frac{f(x_i)}{\prod_{j \ne i}(x_i - x_j)} $$ である。 補間多項式は点の集合だけで決まる（）ので、並べ方によりません。Lagrange 形の $\ell_i$ の $x^k$ の係数は $1/\prod_{j \ne i}(x_i - x_j)$ なので、$\sum_i f(x_i)\ell_i$ の $x^k$ の係数は右辺です。 $x_0, \ldots, x_n$ を相異なる点とする。 1. $f$ を補間する $p_n \in \mathcal{P}_n$ は $$ p_n(x) = \sum_{k=0}^{n} f[x_0, \ldots, x_k]\prod_{j=0}^{k-1}(x - x_j) $$ と書ける（**Newton 形**。$k = 0$ の積は $1$）。 2. $k \ge 1$ について $$ f[x_0, \ldots, x_k] = \frac{f[x_1, \ldots, x_k] - f[x_0, \ldots, x_{k-1}]}{x_k - x_0} $$ が成り立つ。 **1.** $p_k = p_{x_0, \ldots, x_k}$ とおき、$p_k = p_{k-1} + f[x_0, \ldots, x_k]\prod_{j \lt k}(x - x_j)$（$k \ge 1$）を示せば、$k$ について足し合わせて主張を得ます。$r = p_k - p_{k-1}$ は $\mathcal{P}_k$ の元で $x_0, \ldots, x_{k-1}$ で $0$ になり、$p_{k-1}$ の次数は $k-1$ 以下なので $x^k$ の係数は $f[x_0, \ldots, x_k]$ です。$s = r - f[x_0, \ldots, x_k]\prod_{j \lt k}(x - x_j)$ は $x^k$ の項が打ち消し合うので $\mathcal{P}_{k-1}$ の元で、$k$ 個の点 $x_0, \ldots, x_{k-1}$ で $0$ になります。 から $s = 0$ です。 **2.** $q_0 = p_{x_0, \ldots, x_{k-1}}$、$q_1 = p_{x_1, \ldots, x_k}$（ともに $\mathcal{P}_{k-1}$ の元）とし、 $$ P(x) = \frac{(x - x_0)\,q_1(x) - (x - x_k)\,q_0(x)}{x_k - x_0} $$ とおきます。$P \in \mathcal{P}_k$ で、$P(x_0) = q_0(x_0) = f(x_0)$、$P(x_k) = q_1(x_k) = f(x_k)$、$0 \lt i \lt k$ では $P(x_i) = \bigl((x_i - x_0) - (x_i - x_k)\bigr)f(x_i)/(x_k - x_0) = f(x_i)$ です。一意性（）から $P = p_k$ で、両辺の $x^k$ の係数を比べると漸化式を得ます。 漸化式 2 から、差分商は表の形で次々に計算できます（$f[x_i, x_{i+1}]$ の列、$f[x_i, x_{i+1}, x_{i+2}]$ の列、…）。$k = 1, 2$ の場合は、[非線形方程式の反復解法](/mathematics/numerical-analysis/nonlinear-equations)で割線法の解析に使った差分商と一致します。そこで示した 2 階の場合（）は、 の $n = 2$ の場合です。 $I$ を開区間、$k \ge 1$ とし、$g : I \to \mathbb{R}$ は $k$ 回微分可能で、相異なる $k + 1$ 個の点 $t_0 \lt t_1 \lt \cdots \lt t_k$（$I$ の点）で $0$ になるとする。このとき $g^{(k)}(\xi) = 0$ を満たす $\xi \in (t_0, t_k)$ が存在する。 $k$ についての帰納法です。$g$ は微分可能なので各 $[t_{i}, t_{i+1}]$ で連続かつ内部で微分可能で、両端で $0$ なので、Rolle の定理（）から $g'(s_i) = 0$ となる $s_i \in (t_i, t_{i+1})$ があります（$i = 0, \ldots, k-1$）。$s_0 \lt \cdots \lt s_{k-1}$ は相異なり $(t_0, t_k)$ に入ります。$k = 1$ ならこれで終わりです。$k \ge 2$ なら $g'$ は $k-1$ 回微分可能で $k$ 個の点で $0$ なので、帰納法の仮定から $(g')^{(k-1)} = g^{(k)}$ は $(s_0, s_{k-1}) \subset (t_0, t_k)$ の点で $0$ になります。 $n \ge 1$、$I$ を開区間、$f : I \to \mathbb{R}$ を $n$ 回微分可能とし、$x_0, \ldots, x_n \in I$ を相異なる点とする。このとき $\min_i x_i \lt \xi \lt \max_i x_i$ を満たす $\xi$ があって $$ f[x_0, \ldots, x_n] = \frac{f^{(n)}(\xi)}{n!} $$ が成り立つ。 $g = f - p_n$ は $I$ 上 $n$ 回微分可能で、$n+1$ 個の点 $x_0, \ldots, x_n$ で $0$ です。 から $g^{(n)}(\xi) = 0$ となる $\xi$ が最小の点と最大の点の間にあります。$p_n$ の $x^n$ の係数は $f[x_0, \ldots, x_n]$ なので $p_n^{(n)} = n!\,f[x_0, \ldots, x_n]$ で、$g^{(n)}(\xi) = 0$ が主張です。 ## 4. 補間誤差 $I$ を開区間、$f : I \to \mathbb{R}$ を $n+1$ 回微分可能とし、$x_0, \ldots, x_n \in I$ を相異なる点、$p_n \in \mathcal{P}_n$ をそこでの $f$ の補間多項式とする。任意の $x \in I$ について、$x, x_0, \ldots, x_n$ の最小値と最大値の間の $\xi$ があって $$ f(x) - p_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\,\omega(x), \qquad \omega(x) = \prod_{j=0}^{n}(x - x_j) $$ が成り立つ。 証明の方針は、誤差を $\omega(x)$ の定数倍で打ち消す補助関数 $g$ を作り、$g$ の零点が多いことから Rolle の定理（）で $g^{(n+1)}$ の零点を取り出す、というものです。なお「モニック」は最高次の係数が $1$ の多項式のことで、$\omega$ はその例です。 $x$ がどれかの $x_j$ に等しければ両辺は $0$ です。そうでないとき $\omega(x) \ne 0$ なので、$c = (f(x) - p_n(x))/\omega(x)$ とおき、$g(t) = f(t) - p_n(t) - c\,\omega(t)$ を考えます。$g$ は $n+1$ 回微分可能で、$n+2$ 個の相異なる点 $x_0, \ldots, x_n, x$ で $0$ です。（$k = n+1$）から、それらの最小値と最大値の間の $\xi$ で $g^{(n+1)}(\xi) = 0$ です。$p_n^{(n+1)} = 0$、$\omega$ はモニックな $n+1$ 次式なので $\omega^{(n+1)} = (n+1)!$ で、$0 = f^{(n+1)}(\xi) - c(n+1)!$ です。$c$ の定義に戻せば主張を得ます。 誤差は 2 つの因子の積です。$f^{(n+1)}/(n+1)!$ は関数の性質で、こちらは選べません。$\omega(x)$ は点の配置だけで決まり、こちらは選べます。区間の上で $\max|\omega|$ を最小にする点の選び方が、後の節の Chebyshev 点です（）。 $x_0, \ldots, x_n$ を相異なる実数とする。 1. 任意の実数 $y_i, z_i$（$i = 0, \ldots, n$）について、$H(x_i) = y_i$、$H'(x_i) = z_i$ を満たす $H \in \mathcal{P}_{2n+1}$ がただ 1 つ存在する。 2. $I$ を $x_0, \ldots, x_n$ を含む開区間、$f : I \to \mathbb{R}$ を $2n+2$ 回微分可能とし、$H$ を $y_i = f(x_i)$、$z_i = f'(x_i)$ についての 1 の多項式とする。任意の $x \in I$ について、$x, x_0, \ldots, x_n$ の最小値と最大値の間の $\xi$ があって $$ f(x) - H(x) = \frac{f^{(2n+2)}(\xi)}{(2n+2)!}\,\omega(x)^2 $$ が成り立つ。 証明の方針は、1 では未知数と条件の数がどちらも $2n+2$ なので、すべて $0$ のデータから $H = 0$ しか出ないことを示せば十分だ、というものです。2 は補間誤差の公式と同じ手順で、各 $x_i$ が値と傾きの 2 回分の零点として数えられるため、零点が $2n+2$ 個になって $\omega(x)^2$ が現れます。 **1.** $H(x) = \sum_{k=0}^{2n+1} a_k x^k$ とおくと、$2n+2$ 個の条件は係数 $a_0, \ldots, a_{2n+1}$ についての $2n+2$ 元の連立 1 次方程式です。正則性の判定（ の 4 ⇒ 5）から、斉次方程式（$y_i = z_i = 0$）の解が $H = 0$ だけであることを示せば十分です。$H \in \mathcal{P}_{2n+1}$ が $H(x_i) = H'(x_i) = 0$ を満たすとします。点を小さい順に並べると、Rolle の定理（）から隣り合う 2 点の間に $H'$ の零点が 1 つずつ、計 $n$ 個あり、これらは $x_i$ と異なります。$H'(x_i) = 0$ と合わせて、$H' \in \mathcal{P}_{2n}$ は相異なる $2n+1$ 個の零点をもつので、 から $H' = 0$ です。よって $H$ は定数で、$H(x_0) = 0$ から $H = 0$ です。 **2.** $x$ がどれかの $x_j$ なら両辺は $0$ です。そうでないとき $c = (f(x) - H(x))/\omega(x)^2$ とおき、$g(t) = f(t) - H(t) - c\,\omega(t)^2$ を考えます。$g$ は $n + 2$ 個の相異なる点 $x_0, \ldots, x_n, x$ で $0$ なので、Rolle の定理から、これらを小さい順に並べたときの隣り合う 2 点の間に $g'$ の零点が 1 つずつ、計 $n+1$ 個あります。さらに $\omega(x_i) = 0$ から $(\omega^2)'(x_i) = 2\omega(x_i)\omega'(x_i) = 0$ なので $g'(x_i) = f'(x_i) - H'(x_i) = 0$ です。$x_i$ は並べた点そのものなので、点の間にある先の $n+1$ 個とは異なります。よって $g'$ は相異なる $2n+2$ 個の零点をもち、（$g'$ に $k = 2n+1$ で適用）から $g^{(2n+2)}(\xi) = 0$ となる $\xi$ が、$x, x_0, \ldots, x_n$ の最小値と最大値の間にあります。$H^{(2n+2)} = 0$、$\omega^2$ はモニックな $2n+2$ 次式なので $(\omega^2)^{(2n+2)} = (2n+2)!$ で、$f^{(2n+2)}(\xi) = c\,(2n+2)!$ です。 Hermite 補間は、$n+1$ 点の値と傾きという $2n+2$ 個の情報から $2n+1$ 次式を決めます。誤差公式の $\omega(x)^2$ は区間の上で符号を変えないので、積分の平均値の定理と組み合わせると「$f^{(2n+2)}(\xi)/(2n+2)!$ × $\int\omega^2$」の形の誤差公式が得られます。これが [Newton–Cotes 則と Gauss 求積](/mathematics/numerical-analysis/newton-cotes-and-gauss-quadrature)の章で Gauss 求積の誤差を導く道筋です。 $f(x) = 1/(1+25x^2)$ を $[-1,1]$ の $n+1$ 個の等間隔点と、$T_{n+1}$ の零点 $\cos\bigl((2k+1)\pi/(2n+2)\bigr)$（後の節で定義する Chebyshev 点）で補間し、4001 点の格子上の最大誤差を比べます。補間多項式の係数は求めず、SciPy の重心形の補間（）で値を計算します。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 import numpy as np from numpy.polynomial import chebyshev as C from scipy.interpolate import BarycentricInterpolator as BI f = lambda x: 1 / (1 + 25 * x**2) # Runge の関数 xx = np.linspace(-1, 1, 4001) # 誤差を測る点 def err(nodes, g=f): return np.max(np.abs(BI(nodes, g(nodes))(xx) - g(xx))) print(" n 等間隔点 Chebyshev点 chebinterpolate x^5再現") for n in [5, 10, 20, 40]: eq = np.linspace(-1, 1, n + 1) k = np.arange(n + 1) ch = np.cos((2 * k + 1) * np.pi / (2 * n + 2)) # T_{n+1} の零点 c = C.chebinterpolate(f, n) # 同じ点で別の算法 e2 = np.max(np.abs(C.chebval(xx, c) - f(xx))) ok = err(ch, lambda x: x**5) < 1e-12 # 5 次式は誤差 0 のはず print(f"{n:3d} {err(eq):.2e} {err(ch):.2e} {e2:.2e}" f" {ok}") ``` ```text n 等間隔点 Chebyshev点 chebinterpolate x^5再現 5 4.33e-01 5.56e-01 5.56e-01 True 10 1.92e+00 1.09e-01 1.09e-01 True 20 5.98e+01 1.53e-02 1.53e-02 True 40 1.05e+05 2.89e-04 2.89e-04 True ``` **検算.** (a) $n \ge 5$ の補間は 5 次式を正確に再現するはずで（ の一意性）、最後の列は誤差が $10^{-12}$ 未満であることを示しています。(b) Chebyshev 点の誤差は、重心形の計算と、NumPy の `chebinterpolate`（同じ点での補間を Chebyshev 級数の係数として求める別の算法）とで表示の桁まで一致します。(c) 等間隔点の誤差は $n$ を 2 倍にするごとに大きくなり（Runge 現象）、Chebyshev 点の誤差は $n$ を 10 増やすごとにほぼ一定の比で減ります。この比が $\rho^{-10}$（$\rho = (1+\sqrt{26})/5$）であることは、後の例（）で確かめます。 $n = 5$ では等間隔点のほうがわずかに良いことにも注意してください。Chebyshev 点が有利なのは点数を増やしたときの振る舞いで、どの $n$ でも勝つわけではありません。等間隔点で何が起きているのかを、次の節で補間という操作そのものの性質として捉えます。 ## 5. Lebesgue 定数 この節では、補間という操作が入力のずれをどれだけ増幅するかを 1 つの数で測ります。$[a, b]$ 上の連続関数の空間 $C[a,b]$ に[一様ノルム](/mathematics/real-analysis/uniform-convergence#def-sup-norm) $\|g\| = \max_{[a,b]}|g|$（関数の絶対値の最大値）を入れます。$[a,b]$ の相異なる点 $x_0, \ldots, x_n$ を固定し、$f \in C[a,b]$ にその補間多項式を対応させる写像を $\mathcal{I}_n f = \sum_i f(x_i)\ell_i$ と書きます。$\mathcal{I}_n$ は線形で、 の一意性から $q \in \mathcal{P}_n$ について $\mathcal{I}_n q = q$ です。 $\lambda_n(x) = \sum_{i=0}^{n}|\ell_i(x)|$ を点 $x_0, \ldots, x_n$ の **Lebesgue 関数**、$\Lambda_n = \max_{x \in [a,b]}\lambda_n(x)$ を **Lebesgue 定数**という。また $f \in C[a,b]$ について $E_n(f) = \inf_{q \in \mathcal{P}_n}\|f - q\|$ を、$n$ 次以下の多項式による**最良一様近似の誤差**という。 $\lambda_n$ は連続なので、最大値は存在します（）。$\sum_i \ell_i = 1$（定数 $1$ の補間）なので $\lambda_n \ge 1$、したがって $\Lambda_n \ge 1$ です。 次の定理の 1 の「作用素ノルム」は、$\|f\| = 1$ の関数 $f$ に補間を施したとき、結果の大きさが最大でどこまで大きくなるかを表す数のことです。1 はそれが $\Lambda_n$ に等しいと述べていて、証明の後半は、その最大を実現する折れ線の関数を実際に作る部分です。 $[a,b]$ の相異なる点 $x_0, \ldots, x_n$ と $f \in C[a,b]$ について次が成り立つ。 1. $\|\mathcal{I}_n f\| \le \Lambda_n\|f\|$ であり、等号を満たす $f \ne 0$ が存在する。すなわち $\Lambda_n$ は $\mathcal{I}_n$ の作用素ノルムである。 2. $\|f - \mathcal{I}_n f\| \le (1 + \Lambda_n)\,E_n(f)$。 3. $E_n(f) \to 0$（$n \to \infty$）。 **1.** $|\mathcal{I}_n f(x)| \le \sum_i |f(x_i)||\ell_i(x)| \le \|f\|\lambda_n(x) \le \|f\|\Lambda_n$ です。等号のために、$\lambda_n(\hat x) = \Lambda_n$ となる $\hat x$ を取り、$s_i$ を $\ell_i(\hat x)$ の符号（$0$ なら $1$）とします。点を小さい順に並べ、$f$ を各点で $s_i$ をとる折れ線（両端の点の外では一定）とすると、$f$ は連続で $\|f\| = 1$ であり、$\mathcal{I}_n f(\hat x) = \sum_i s_i\ell_i(\hat x) = \lambda_n(\hat x) = \Lambda_n$ です。 **2.** 任意の $q \in \mathcal{P}_n$ について $\mathcal{I}_n q = q$ なので、$f - \mathcal{I}_n f = (f - q) - \mathcal{I}_n(f - q)$ です。1 から $\|f - \mathcal{I}_n f\| \le (1 + \Lambda_n)\|f - q\|$ で、$q$ について下限をとれば主張を得ます。 **3.** $\varepsilon \gt 0$ について、Weierstrass の多項式近似定理（）から $\|f - P\| \lt \varepsilon$ となる多項式 $P$ があります。$P$ の次数を $N$ とすると、$n \ge N$ で $P \in \mathcal{P}_n$ なので $E_n(f) \le \varepsilon$ です。 2 と 3 を合わせても、$\|f - \mathcal{I}_n f\| \to 0$ とは言えません。$E_n(f)$ が $0$ に近づく速さより $\Lambda_n$ が増える速さのほうが速ければ、補間は発散しえます。Runge の現象はまさにこの状況です。 データ $f(x_i)$ に大きさ $\varepsilon$ 以下の誤差（測定誤差や丸め誤差）$\delta_i$ が乗ると、補間多項式は $\sum_i\delta_i\ell_i$ だけずれます。その一様ノルムは $\varepsilon\Lambda_n$ 以下で、 の 1 から、最悪の場合（$\delta_i = \varepsilon s_i$）にはちょうど $\varepsilon\Lambda_n$ になります。つまり $\Lambda_n$ は「データから補間多項式への写像」の（一様ノルムで測った絶対）条件数です。倍精度の丸め誤差 $u \approx 1.1 \times 10^{-16}$ を想定すると、$\Lambda_n \approx 10^{10}$ の点の配置では補間の値は最悪で 6 桁程度しか信用できません。後の表（）では等間隔点の $\Lambda_{30}$ が約 $6.6 \times 10^6$ で、$n$ を 1 増やすごとにおよそ 2 倍になっていくので、$n = 40$ では約 $5 \times 10^9$ とこの水準に迫ります。 $n \ge 1$、$h = (b-a)/n$、$x_j = a + jh$（$j = 0, \ldots, n$）とする。このとき $$ \Lambda_n \ge \lambda_n\Bigl(a + \frac h2\Bigr) \ge \frac{2^n}{4n^2} $$ である。 証明の方針は、端の最初の小区間の中点 $a + h/2$ だけに注目し、そこでの各 $|\ell_j|$ を下から評価して足し合わせる、というものです。足し合わせには二項定理を使うので、式変形が続きますが、最後は $\binom{2n}{n}$ の下からの評価に帰着します。 $x = a + th$ とおくと $\ell_j(x) = \prod_{k \ne j}(t - k)/(j - k)$ です。$t = 1/2$ では $$ \prod_{k=0}^{n}\Bigl|\frac12 - k\Bigr| = \frac12\prod_{k=1}^{n}\frac{2k - 1}{2} = \frac{(2n-1)!!}{2^{n+1}}, \qquad \prod_{k \ne j}|j - k| = j!\,(n-j)! $$ なので（$(2n-1)!! = 1 \cdot 3 \cdots (2n-1)$）、 $$ \Bigl|\ell_j\Bigl(a + \frac h2\Bigr)\Bigr| = \frac{(2n-1)!!}{2^{n+1}\,|1/2 - j|\; j!\,(n-j)!} \ge \frac{(2n-1)!!}{2^{n+1}\, n\; j!\,(n-j)!} $$ です（$|1/2 - j| \le n$ を使いました）。二項定理から $\sum_{j=0}^{n} 1/(j!(n-j)!) = 2^n/n!$ なので $$ \lambda_n\Bigl(a + \frac h2\Bigr) \ge \frac{(2n-1)!!}{2n \cdot n!} = \frac{(2n)!}{2^{n+1}\,n\,(n!)^2} = \frac{1}{2^{n+1}n}\binom{2n}{n} $$ です（$(2n-1)!! = (2n)!/(2^n n!)$）。最後に $\binom{2n}{n} \ge 4^n/(2n)$ を示します。$\binom{2n}{k+1}/\binom{2n}{k} = (2n-k)/(k+1)$ は $k \le n-1$ で $1$ 以上、$k \ge n$ で $1$ 未満なので、$\binom{2n}{n}$ は $\binom{2n}{k}$ の最大です。二項定理から $4^n = \sum_{k=0}^{2n}\binom{2n}{k} = 2 + \sum_{k=1}^{2n-1}\binom{2n}{k} \le 2 + (2n-1)\binom{2n}{n} \le 2n\binom{2n}{n}$（$\binom{2n}{n} \ge 2$ を使いました）です。合わせて $\lambda_n(a + h/2) \ge 4^n/(2^{n+1}n \cdot 2n) = 2^n/(4n^2)$ です。 下界は $n = 10$ で $2.56$、$n = 30$ で約 $3.0 \times 10^5$ と控えめですが、指数的に増えることが要点です。実際の値はもっと大きく（後の表で $n = 30$ で約 $6.6 \times 10^6$）、およそ $2^{n+1}/(e\,n\log n)$ の速さで増えることが知られています（Trefethen『Approximation Theory and Approximation Practice』第 15 章）。 ## 6. Chebyshev 多項式と Chebyshev 点 Chebyshev 多項式 $T_k$ は、漸化式 $T_0 = 1$、$T_1 = x$、$T_{k+1} = 2xT_k - T_{k-1}$ で定めます（）。 この多項式の要点は、$x = \cos\theta$ とおくと $T_k(x) = \cos k\theta$ になることです。角度 $\theta$ を等間隔に動かすと $\cos k\theta$ は等しい振れ幅で振動し、$[-1,1]$ 上で $|T_k| \le 1$ が保たれます。零点も角度が等間隔の点に対応するので、**Chebyshev 点**は、半円周上に等間隔に並べた点を横軸へ落とした点になります。下の図は $n = 5$（6 点）の場合です。
-1 0 1
$k \ge 0$ について次が成り立つ。 1. すべての実数 $\theta$ について $T_k(\cos\theta) = \cos k\theta$。とくに $[-1,1]$ 上で $|T_k| \le 1$。 2. $k \ge 1$ なら $T_k$ は $k$ 次で、最高次の係数は $2^{k-1}$ である。$T_k(-x) = (-1)^kT_k(x)$。 3. $k \ge 1$ なら $T_k$ の零点は $\cos\bigl((2j+1)\pi/(2k)\bigr)$（$j = 0, \ldots, k-1$）の $k$ 個で、すべて $(-1, 1)$ の相異なる点である。また $\eta_j = \cos(j\pi/k)$（$j = 0, \ldots, k$）で $T_k(\eta_j) = (-1)^j$ である。 **1.** $k = 0, 1$ では定義どおりです。$k - 1, k$ で成り立てば、加法定理 $\cos(k+1)\theta + \cos(k-1)\theta = 2\cos\theta\cos k\theta$ から $$ T_{k+1}(\cos\theta) = 2\cos\theta\,T_k(\cos\theta) - T_{k-1}(\cos\theta) = 2\cos\theta\cos k\theta - \cos(k-1)\theta = \cos(k+1)\theta $$ です。$x \in [-1,1]$ は $x = \cos\theta$ と書けるので $|T_k(x)| \le 1$ です。 **2.** 漸化式から帰納法で、$T_{k+1}$ の最高次の項は $2x \cdot 2^{k-1}x^k = 2^k x^{k+1}$ です（$k = 0$ から $1$ へは $T_1 = x$）。偶奇も漸化式から帰納法で従います。 **3.** $\theta_j = (2j+1)\pi/(2k) \in (0, \pi)$ は相異なり、$\cos$ は $[0, \pi]$ で狭義減少なので $\cos\theta_j$ は $(-1,1)$ の相異なる $k$ 点です。1 から $T_k(\cos\theta_j) = \cos\bigl((2j+1)\pi/2\bigr) = 0$ です。ほかにも零点があれば、$T_k \in \mathcal{P}_k$ が $k+1$ 個の相異なる零点をもつので から $T_k = 0$ となり、最高次の係数が $2^{k-1} \ne 0$ であることに矛盾します。よってこれがすべての零点です。$T_k(\cos(j\pi/k)) = \cos j\pi = (-1)^j$ も 1 から従います。 $n \ge 0$ とする。 1. $[-1,1]$ 上で、$n+1$ 次のモニック多項式 $q$ は $\max_{[-1,1]}|q| \ge 2^{-n}$ を満たし、$q = 2^{-n}T_{n+1}$ で等号が成り立つ。 2. $a \lt b$ とする。$[a,b]$ 上で、$n+1$ 次のモニック多項式 $q$ は $\max_{[a,b]}|q| \ge 2\bigl((b-a)/4\bigr)^{n+1}$ を満たす。 **1.** $\tilde T = 2^{-n}T_{n+1}$ は の 2 からモニックで、1 から $\max|\tilde T| \le 2^{-n}$、3 から $\eta_j = \cos(j\pi/(n+1))$（$j = 0, \ldots, n+1$）で $\tilde T(\eta_j) = (-1)^j2^{-n}$ です。よって $\max|\tilde T| = 2^{-n}$ です。 $\max|q| \lt 2^{-n}$ を満たすモニックな $q$ があったと仮定します。$r = \tilde T - q$ は最高次の項が打ち消し合うので $\mathcal{P}_n$ の元です。$|q(\eta_j)| \lt 2^{-n}$ なので $r(\eta_j) = (-1)^j2^{-n} - q(\eta_j)$ の符号は $(-1)^j$ です。$\eta_0 \gt \eta_1 \gt \cdots \gt \eta_{n+1}$ の隣り合う 2 点で $r$ の符号は逆なので、中間値の定理（）から各区間 $(\eta_{j+1}, \eta_j)$ に $r$ の零点があり、相異なる零点が $n+1$ 個得られます。 から $r = 0$ ですが、$r(\eta_0) \gt 0$ に矛盾します。 **2.** $x = \frac{a+b}{2} + \frac{b-a}{2}t$ とおくと、$Q(t) = (2/(b-a))^{n+1}q\bigl(\frac{a+b}{2} + \frac{b-a}{2}t\bigr)$ は $t$ の $n+1$ 次のモニック多項式です。1 から $\max_{t \in [-1,1]}|Q(t)| \ge 2^{-n}$ なので、$\max_{[a,b]}|q| \ge ((b-a)/2)^{n+1}2^{-n} = 2((b-a)/4)^{n+1}$ です。 この定理を補間誤差の公式に組み合わせると、$\max|\omega|$ を最小にする点が決まります。 $a \lt b$、$n \ge 0$ とし、**Chebyshev 点** $$ x_j = \frac{a+b}{2} + \frac{b-a}{2}\cos\frac{(2j+1)\pi}{2n+2} \qquad (j = 0, \ldots, n) $$ をとる。このとき $\max_{[a,b]}|\omega| = 2((b-a)/4)^{n+1}$ で、これは相異なる $n+1$ 点のすべての選び方の中で最小である。$f$ が $[a,b]$ を含む開区間で $n+1$ 回微分可能で $|f^{(n+1)}| \le M_{n+1}$（$[a,b]$ 上）ならば $$ \max_{[a,b]}|f - p_n| \le \frac{2M_{n+1}}{(n+1)!}\Bigl(\frac{b-a}{4}\Bigr)^{n+1} $$ である。とくに $[a,b] = [-1,1]$ では右辺は $M_{n+1}/(2^n(n+1)!)$ である。 $t(x) = (2x - a - b)/(b-a)$ とおきます。 の 2・3 から、$x$ の多項式 $T_{n+1}(t(x))$ は $n+1$ 次で、最高次の係数は $2^n(2/(b-a))^{n+1}$、零点はちょうど $x_0, \ldots, x_n$ です。$\omega - T_{n+1}(t(x))/\bigl(2^n(2/(b-a))^{n+1}\bigr)$ は $\mathcal{P}_n$ の元で $n+1$ 個の点 $x_j$ で $0$ なので、 から $\omega(x) = 2^{-n}((b-a)/2)^{n+1}T_{n+1}(t(x))$ です。$t(x)$ は $[a,b]$ を $[-1,1]$ に写すので、 の 1・3 から $\max_{[a,b]}|\omega| = 2^{-n}((b-a)/2)^{n+1} = 2((b-a)/4)^{n+1}$ です。他の点の選び方でも $\omega$ はモニックな $n+1$ 次式なので、 の 2 からこれ以上です。誤差の評価は、 の $\xi$ が $[a,b]$ に入ることから従います。 $x \ge 1$、$w = x + \sqrt{x^2 - 1}$ とすると、$k \ge 0$ について $$ T_k(x) = \frac{w^k + w^{-k}}{2}, \qquad \frac{w^k}{2} \le T_k(x) \le w^k $$ である。また $x \le -1$ では $|T_k(x)| = T_k(|x|)$ である。 $w \ge 1$ で $w^{-1} = x - \sqrt{x^2 - 1}$ なので $w + w^{-1} = 2x$ です。$S_k = (w^k + w^{-k})/2$ とおくと $S_0 = 1$、$S_1 = x$ で、 $$ 2xS_k - S_{k-1} = \frac{(w + w^{-1})(w^k + w^{-k}) - w^{k-1} - w^{1-k}}{2} = \frac{w^{k+1} + w^{-k-1}}{2} = S_{k+1} $$ なので、$T_k$ と同じ漸化式と初期値を満たし、帰納法で $T_k(x) = S_k$ です。$0 \lt w^{-k} \le 1 \le w^k$ から不等式が従います。最後の主張は の 2 の偶奇性です。 この補題は、区間 $[-1,1]$ の中では $1$ 以下に抑えられている $T_k$ が、外に出ると $w^k$ の速さで増えることを述べています。「区間の上では小さく、区間の外の 1 点では大きい」多項式が必要な場面（[共役勾配法](/mathematics/numerical-linear-algebra/lanczos-and-conjugate-gradient#def-cg)の収束の評価など）で、Chebyshev 多項式が最適な道具になる理由です（ を参照）。 $n \ge 0$ とし、$[-1,1]$ の Chebyshev 点 $x_j = \cos\bigl((2j+1)\pi/(2n+2)\bigr)$（$j = 0, \ldots, n$）の Lebesgue 定数を $\Lambda_n$ とする。このとき $$ \Lambda_n \le \frac{2}{\pi}\log(2n+1) + 1 + \frac{4}{\pi} $$ である。$[a,b]$ の Chebyshev 点（）でも同じ評価が成り立つ。 証明は長いので Appendix に置きます。最後の主張は、$\ell_j$ が 1 次の変数変換 $t(x)$ で $[-1,1]$ の場合の $\ell_j$ に移ることから従います。等間隔点の指数的な増大と比べると、Chebyshev 点の Lebesgue 定数は $n = 1000$ でもこの評価から $7.2$ 以下です。 格子（各小区間を 20 等分した点）の上で Lebesgue 関数の最大を求め、等間隔点では の下界と、Chebyshev 点では の上界と比べます。等間隔点の値は $\ell_j$ の積に大きな数と小さな数が混ざるので、丸めの影響でないことを mpmath の 40 桁の計算で確かめます。格子上の最大は $\Lambda_n$ を下から近似する値です。 ```python # Python 3.12.10 / NumPy 2.5.2 / mpmath 1.3.0 import numpy as np from mpmath import mp, mpf def leb_np(nd, xs): # 格子 xs 上の Lebesgue 関数の最大 L = np.ones((len(nd), len(xs))) for j in range(len(nd)): for k in range(len(nd)): if k != j: L[j] *= (xs - nd[k]) / (nd[j] - nd[k]) return np.abs(L).sum(axis=0).max() def leb_mp(nd, xs): # 同じ量を mpmath の 40 桁で mp.dps = 40 nd, best = [mpf(t) for t in nd], mpf(0) for x in map(mpf, xs): s = sum(abs(mp.fprod((x - b) / (a - b) for b in nd if b != a)) for a in nd) best = max(best, s) return best print(" n 等間隔(倍) 等間隔(40桁) 下界 Cheb 上界") for n in [5, 10, 20, 30]: xs = np.linspace(-1, 1, 20 * n + 1) # 各小区間を 20 等分 eq = np.linspace(-1, 1, n + 1) ch = np.cos((2 * np.arange(n + 1) + 1) * np.pi / (2 * n + 2)) lo = 2.0**n / (4 * n * n) up = 2 / np.pi * np.log(2 * n + 1) + 1 + 4 / np.pi print(f"{n:3d} {leb_np(eq, xs):.4e} {float(leb_mp(eq, xs)):.4e}" f" {lo:.2e} {leb_np(ch, xs):.3f} {up:.3f}") ``` ```text n 等間隔(倍) 等間隔(40桁) 下界 Cheb 上界 5 3.1050e+00 3.1050e+00 3.20e-01 2.104 3.800 10 2.9891e+01 2.9891e+01 2.56e+00 2.489 4.211 20 1.0987e+04 1.0987e+04 6.55e+02 2.901 4.637 30 6.5710e+06 6.5710e+06 2.98e+05 3.149 4.890 ``` **検算.** (a) $n = 2$ の場合は手で計算でき（）、等間隔点で $5/4$、Chebyshev 点で $5/3$ です。(b) 等間隔点の値は倍精度と 40 桁で表示の 5 桁まで一致し、丸め誤差の影響ではありません。(c) 等間隔点の値は下界 $2^n/(4n^2)$ の上にあり、$n$ を 10 増やすごとに数百倍（$20 \to 30$ で約 $600$ 倍、1 増やすごとにおよそ $1.9$ 倍）になって、下界の $2^n$ に近い速さで増えています。Chebyshev 点の値は上界の下にあり、$n$ を 10 増やしても $0.25$〜$0.41$ ずつしか増えません。参考までに、よく引用される評価 $(2/\pi)\log(n+1) + 1$ は $n = 30$ で $3.186$ で、表の $3.149$ はその下にもあります。 Lebesgue 定数と最良近似を組み合わせると、Chebyshev 点での補間が解析関数に対して幾何級数的に収束することが、導関数を使わずに示せます。 次の系の「Bernstein 楕円 $E_\rho$」は、複素数の平面で $[-1,1]$ を囲む、焦点が $\pm1$ の楕円で、$\rho$ が大きいほど大きな楕円です。「正則」は複素数の変数で微分できることで、$f$ が $E_\rho$ の内側でこの意味でなめらかなら、$\rho$ が大きいほど補間の誤差が速く減る、というのが系の内容です。 $\rho \gt 1$ とし、$f$ は Bernstein 楕円 $E_\rho$ で正則、$[-1,1]$ 上で実数値とする。$p_n$ を $[-1,1]$ の Chebyshev 点での $f$ の補間多項式とすると、$1 \lt r \lt \rho$ と $M_r = \max_{\partial E_r}|f|$ について $$ \max_{[-1,1]}|f - p_n| \le \Bigl(2 + \frac4\pi + \frac2\pi\log(2n+1)\Bigr)\frac{2M_r\,r^{-n}}{r - 1} $$ が成り立つ。 Chebyshev 展開の打ち切り（。係数の減衰は ）$\sum_{k=0}^{n}a_kT_k$ は $\mathcal{P}_n$ の元で、$f$ との差は $[-1,1]$ 上で $2M_rr^{-n}/(r-1)$ 以下です。よって $E_n(f) \le 2M_rr^{-n}/(r-1)$ です。 の 2 と から主張を得ます。 対数の因子は幾何級数の減衰に比べて無視できるので、誤差はおよそ $r^{-n}$ の速さで減ります。$r$ は $\rho$ にいくらでも近く取れるので、収束の速さを決めるのは $f$ が正則でいられる最大の Bernstein 楕円の大きさ $\rho$ です。Runge の関数は $\pm i/5$ に極をもち、$\rho = 1/5 + \sqrt{1 + 1/25} = (1 + \sqrt{26})/5 \approx 1.2198$ です。 $e^x$ には の評価 $e/(2^n(n+1)!)$ を、Runge の関数には の予言する速さ $\rho^{-n}$ を当てはめます。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 import math import numpy as np from scipy.interpolate import BarycentricInterpolator as BI xx = np.linspace(-1, 1, 4001) def cheb_err(g, n): # T_{n+1} の零点で補間 k = np.arange(n + 1) x = np.cos((2 * k + 1) * np.pi / (2 * n + 2)) return np.max(np.abs(BI(x, g(x))(xx) - g(xx))) print(" n e^x の誤差 評価 e/(2^n (n+1)!) 比") for n in [2, 4, 6, 8, 10]: e = cheb_err(np.exp, n) b = math.e / (2**n * math.factorial(n + 1)) print(f"{n:3d} {e:.2e} {b:.2e} {e / b:.3f}") rho = (1 + math.sqrt(26)) / 5 # Bernstein 楕円の大きさ runge = lambda x: 1 / (1 + 25 * x**2) print(" n Runge の誤差 誤差×ρ^n") for n in [10, 20, 40, 80, 120]: e = cheb_err(runge, n) print(f"{n:3d} {e:.2e} {e * rho**n:.3f}") print(f"ρ = {rho:.4f}, ρ^(-10) = {rho**-10:.4f}") ``` ```text n e^x の誤差 評価 e/(2^n (n+1)!) 比 2 5.65e-02 1.13e-01 0.499 4 6.40e-04 1.42e-03 0.452 6 3.62e-06 8.43e-06 0.430 8 1.22e-08 2.93e-08 0.417 10 2.71e-11 6.65e-11 0.408 n Runge の誤差 誤差×ρ^n 10 1.09e-01 0.796 20 1.53e-02 0.816 40 2.89e-04 0.819 80 1.02e-07 0.818 120 3.62e-11 0.819 ρ = 1.2198, ρ^(-10) = 0.1371 ``` **検算.** (a) $e^x$ の誤差はすべて評価 $e/(2^n(n+1)!)$ 以下で（比が 1 未満）、比がほぼ一定なので、誤差は評価式と同じ速さで減っています。比が $1$ より小さいのは、評価が $|f^{(n+1)}|$ の最大値 $e$ を使っているためです。(b) Runge の関数の $n = 10, 20, 40$ の誤差は、別の点の配置と並べた先の例（）の Chebyshev 点の列と一致します。そこでは別の算法（`chebinterpolate`）とも一致していました。(c) 誤差に $\rho^n$ を掛けた値は $0.82$ 付近で一定になり、誤差がちょうど $\rho^{-n}$ の速さで減ること、すなわち $n$ を 10 増やすごとに $\rho^{-10} \approx 0.137$ 倍になることを示しています。これは の予言する速さ（任意の $r \lt \rho$ について $r^{-n}$）と合っています。 Lagrange 形を $w_j = 1/\prod_{k \ne j}(x_j - x_k)$ を使って変形すると、$p_n(x) = \sum_j \frac{w_j f(x_j)}{x - x_j} \Big/ \sum_j \frac{w_j}{x - x_j}$ と書けます（分母は定数 $1$ の補間から出ます）。これを**重心形**といい、重み $w_j$ を一度計算すれば 1 点あたり $O(n)$ の演算で値が求まり、点を足すときの更新も簡単です。Chebyshev 点では $w_j$ が閉じた形で書けます。重心形の計算は丸め誤差に対して安定で、誤差は Lebesgue 定数と単位丸めの積程度に抑えられることが N. J. Higham, "The numerical stability of barycentric Lagrange interpolation", *IMA J. Numer. Anal.* 24 (2004), 547–556 で示されています。総説は J.-P. Berrut, L. N. Trefethen, "Barycentric Lagrange interpolation", *SIAM Review* 46 (2004), 501–517（[doi:10.1137/S0036144502417715](https://doi.org/10.1137/S0036144502417715)）です。SciPy の `BarycentricInterpolator` はこの公式を使っています。 $E_n(f)$ の下限は達成され（$\mathcal{P}_n$ は有限次元なので、最良近似多項式 $q^*$ が存在する）、$q^*$ はただ 1 つで、$f - q^*$ が $[a,b]$ の $n+2$ 個の点で符号を交互に変えながら絶対値 $E_n(f)$ をとることで特徴づけられます（Chebyshev の交代定理）。 はこの定理の $f = x^{n+1}$ の場合にあたります。$q^*$ を数値で求める算法が Remez 法です。証明は T. J. Rivlin『An Introduction to the Approximation of Functions』（Dover, 1981）第 1 章、Trefethen『Approximation Theory and Approximation Practice』第 10 章にあります。本章では の 2 のように下限 $E_n(f)$ を直接使うので、存在と一意性は必要ありません。 次数を上げずに精度を上げるもう 1 つの方法は、区間を小区間に分けて各小区間で低次の多項式を使い、つなぎ目で滑らかさを課すことです。3 次スプライン（つなぎ目で $C^2$ 級）は、刻み $h$ について誤差 $O(h^4)$ で、点を増やしても Runge 現象のような不安定さが起こりません。理論と算法は C. de Boor『A Practical Guide to Splines』（Springer, 1978、改訂版 2001）を参照してください。 ## 7. 演習 $f(x) = x^3$ の、点 $0, 1, 2, 3$ での差分商の表をの漸化式で作り、Newton 形が $x^3$ に一致すること、$f[0,1,2,3]$ がの $f'''(\xi)/3!$ に等しいことを確かめてください。 値は $0, 1, 8, 27$、1 階の差分商は $1, 7, 19$、2 階は $(7-1)/2 = 3$、$(19-7)/2 = 6$、3 階は $(6-3)/3 = 1$ です。Newton 形は $$ 0 + 1\cdot x + 3x(x-1) + 1\cdot x(x-1)(x-2) = x + 3x^2 - 3x + x^3 - 3x^2 + 2x = x^3 $$ で一致します。$f''' = 6$ なので $f'''(\xi)/3! = 1 = f[0,1,2,3]$ で、$\xi$ はどこでもかまいません。 $[0,1]$ の Chebyshev 点で $e^x$ を補間します。 の評価で誤差 $10^{-8}$ 以下が保証される最小の $n$ を求めてください。 $[0,1]$ 上 $|f^{(n+1)}| \le e$ なので、評価は $2e(1/4)^{n+1}/(n+1)!$ です。$n = 6$ では $2e/(4^7 \cdot 7!) = 2e/82{,}575{,}360 \approx 6.6 \times 10^{-8}$ で足りず、$n = 7$ では $2e/(4^8\cdot 8!) = 2e/2{,}642{,}411{,}520 \approx 2.1 \times 10^{-9}$ です。よって $n = 7$（8 点）です。 $[-1,1]$ の 3 点について Lebesgue 定数 $\Lambda_2$ を求めてください。(i) 等間隔点 $-1, 0, 1$。(ii) Chebyshev 点 $-\sqrt3/2, 0, \sqrt3/2$。どちらが小さいか、その理由も考えてください。 (i) $\ell_0 = x(x-1)/2$、$\ell_1 = 1 - x^2$、$\ell_2 = x(x+1)/2$ です。対称性から $[0,1]$ で考えると $\lambda_2 = x(1-x)/2 + (1 - x^2) + x(x+1)/2 = 1 + x - x^2$ で、最大は $x = 1/2$ での $5/4$ です。 (ii) $s = \sqrt3/2$ とすると $\ell_0 = x(x-s)/(2s^2)$、$\ell_1 = (s^2 - x^2)/s^2$、$\ell_2 = x(x+s)/(2s^2)$ です。$0 \le x \le s$ では $\lambda_2 = 1 + (xs - x^2)/s^2$ で最大は $x = s/2$ での $5/4$、$s \le x \le 1$ では $\ell_1 \le 0$ なので $\lambda_2 = 2x^2/s^2 - 1$ で、最大は $x = 1$ での $8/3 - 1 = 5/3$ です。よって $\Lambda_2 = 5/3$ です。 $n = 2$ では等間隔点のほうが小さく、Chebyshev 点（第 1 種）は端点 $\pm1$ を含まないので、端の近くで補間が外挿になって大きくなります。Chebyshev 点の利点は $n$ が大きいときの対数的な増え方にあり、数値の例（）では $n = 5$ ですでに逆転しています。 $0 \lt a \lt b$、$\kappa = b/a$ とします。$p(0) = 1$ を満たす $k$ 次以下の多項式 $p$ で $$ \max_{x \in [a,b]}|p(x)| \le 2\Bigl(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\Bigr)^k $$ となるものがあることを、 を使って示してください（この評価は共役勾配法の収束の評価の形です）。 $p(x) = T_k\bigl((b + a - 2x)/(b - a)\bigr)/T_k(x_0)$、$x_0 = (b+a)/(b-a)$ とおきます。$x_0 \gt 1$ なので から $T_k(x_0) \gt 0$ で、$p(0) = 1$ です。$x \in [a,b]$ なら $(b+a-2x)/(b-a) \in [-1,1]$ なので、 の 1 から分子の絶対値は $1$ 以下です。分母について、$x_0 = (\kappa+1)/(\kappa-1)$、$x_0^2 - 1 = 4\kappa/(\kappa-1)^2$ なので $$ w = x_0 + \sqrt{x_0^2 - 1} = \frac{\kappa + 1 + 2\sqrt\kappa}{\kappa - 1} = \frac{(\sqrt\kappa + 1)^2}{(\sqrt\kappa - 1)(\sqrt\kappa + 1)} = \frac{\sqrt\kappa + 1}{\sqrt\kappa - 1} $$ で、 から $T_k(x_0) \ge w^k/2$ です。よって $\max_{[a,b]}|p| \le 2w^{-k}$ です。 ## 参考文献 - L. N. Trefethen, *Approximation Theory and Approximation Practice*, SIAM, 2013（Extended Edition, 2019）— Chapters 3–4（Chebyshev 多項式と補間）、8（解析関数の収束）、13（等間隔点と Runge 現象）、15（Lebesgue 定数）。 - K. E. Atkinson, *An Introduction to Numerical Analysis*, 2nd ed., Wiley, 1989 — Chapter 3（補間理論）、Chapter 4（関数近似）。 - T. J. Rivlin, *Chebyshev Polynomials: From Approximation Theory to Algebra and Number Theory*, 2nd ed., Wiley, 1990 — Chapter 1。 - 山本哲朗『数値解析入門』（増訂版）サイエンス社、2003 — 補間の章。 - J.-P. Berrut, L. N. Trefethen, "Barycentric Lagrange interpolation", *SIAM Review* 46 (2004), 501–517. [doi:10.1137/S0036144502417715](https://doi.org/10.1137/S0036144502417715) ## Appendix: Chebyshev 点の Lebesgue 定数の上界の証明 証明の流れは次のとおりです。まず $\ell_j$ を $\theta$ で表し、次に $\lambda_n$ を余接（$\cot$）の和で上から押さえます。最後に、$x$ に最も近い節点の項だけを別に扱い、残りの項を調和級数に似た和 $\sum 1/(2i-1)$ で評価して、それが対数で抑えられることを使います。 **記号.** $N = n + 1$、$\theta_j = (2j+1)\pi/(2N)$（$j = 0, \ldots, n$）とし、$x \in [-1,1]$ を $x = \cos\theta$（$\theta \in [0,\pi]$）と書きます。証明に使う初等的な不等式を先に確かめます。平均値の定理（）から、$0 \lt u \lt \pi/2$ で $\tan u = u/\cos^2 c \ge u$（$0 \lt c \lt u$）、したがって $\cot u \le 1/u$ です。$\cot$ は $(0, \pi/2]$ で正、導関数 $-1/\sin^2 u \lt 0$ なので狭義減少です（）。$y \gt 0$ で $\log(1+y) = y/(1+c) \ge y/(1+y)$（$0 \lt c \lt y$）です。同じく平均値の定理から、実数 $a, b$ について $|\cos Na - \cos Nb| = N|\sin Nc|\,|a - b| \le N|a - b|$（$c$ は $a$ と $b$ の間）です。 **$\ell_j$ の表示.** の証明のとおり $\omega = 2^{-n}T_N$ で、積の微分から $\omega'(x_j) = \prod_{k \ne j}(x_j - x_k)$ なので、$x \ne x_j$ で $\ell_j(x) = \omega(x)/\bigl((x - x_j)\omega'(x_j)\bigr) = T_N(x)/\bigl((x - x_j)T_N'(x_j)\bigr)$ です。 の 1 の両辺を $\theta$ で微分すると $T_N'(\cos\theta)\sin\theta = N\sin N\theta$ で、$\sin N\theta_j = \sin\bigl((2j+1)\pi/2\bigr) = (-1)^j$ から $|T_N'(x_j)| = N/\sin\theta_j$ です。よって $x = \cos\theta \ne x_j$ で $$ |\ell_j(x)| = \frac{|\cos N\theta|}{N}\cdot\frac{\sin\theta_j}{|\cos\theta - \cos\theta_j|} $$ です。$x$ がどれかの節点なら $\lambda_n(x) = 1$ なので、以下 $x$ は節点でないとします。 **余接への分解.** $\alpha = (\theta_j + \theta)/2$、$\beta = (\theta_j - \theta)/2$ とおくと、$\cos\theta - \cos\theta_j = 2\sin\alpha\sin\beta$、$\sin\theta_j = \sin(\alpha + \beta) = \sin\alpha\cos\beta + \cos\alpha\sin\beta$ なので $$ \frac{\sin\theta_j}{\cos\theta - \cos\theta_j} = \frac12\bigl(\cot\beta + \cot\alpha\bigr) $$ です（$0 \lt \alpha \lt \pi$、$\beta \notin \pi\mathbb{Z}$）。$|\cot|$ は偶関数で周期 $\pi$ なので $|\cot\alpha| = |\cot\bigl((\theta - (2\pi - \theta_j))/2\bigr)|$ です。$\{\theta_j\} \cup \{2\pi - \theta_j\}$ はちょうど $s_m = (2m+1)\pi/(2N)$（$m = 0, \ldots, 2N-1$）の全体なので、$v_m = (\theta - s_m)/2$ とおくと $$ \lambda_n(x) \le \frac{|\cos N\theta|}{2N}\sum_{m=0}^{2N-1}|\cot v_m| $$ です。 **点の並び.** 実数 $v$ と $\pi\mathbb{Z}$ の距離を $\delta(v) \in [0, \pi/2]$ と書くと、$|\cot v| = \cot\delta(v)$ です。$v_m$ は公差 $-h$（$h = \pi/(2N)$、$2Nh = \pi$）の $2N$ 項の等差数列なので、$\pi$ を法として周長 $\pi$ の円周上に間隔 $h$ で並ぶ $2N$ 点です。$\delta(v_{m_0}) = d$ が最小となる $m_0$ を取ると $d \le h/2$ です。$x$ は節点でないので $d \gt 0$ です。$m \ne m_0$ では $v_m - v_{m_0} = kh$ で、$k$ を $2N$ で割った余り $k'$ は $1, \ldots, 2N-1$ を 1 回ずつとり、$\delta(kh) = \min(k', 2N - k')\,h$ です。$v_{m_0} = p\pi + e$（$p \in \mathbb{Z}$、$|e| = d$）と書くと $\delta(v_m) = \delta(e + kh) \ge \delta(kh) - d$ なので、$i = \min(k', 2N - k')$ として $\delta(v_m) \ge (i - 1/2)h$ です。$i$ は $1, \ldots, N-1$ を 2 回ずつ、$N$ を 1 回とります。$\cot$ の減少性から $$ \sum_{m \ne m_0}|\cot v_m| \le 2\sum_{i=1}^{N}\cot\frac{(2i-1)\pi}{4N} \le 2\sum_{i=1}^{N}\frac{4N}{(2i-1)\pi} $$ です。 **最も近い項.** $\cos(Ns_{m_0}) = \cos\bigl((2m_0+1)\pi/2\bigr) = 0$ で、$\theta - (s_{m_0} + 2p\pi) = 2e$ なので、$|\cos N\theta| = |\cos N\theta - \cos N(s_{m_0} + 2p\pi)| \le 2Nd$ です。よって $$ \frac{|\cos N\theta|}{2N}|\cot v_{m_0}| \le \frac{2Nd}{2N}\cot d \le d \cdot \frac1d = 1 $$ です。 **まとめ.** 残りの項は $|\cos N\theta| \le 1$ と上の評価から $\frac{1}{2N}\cdot\frac{8N}{\pi}\sum_{i=1}^{N}\frac{1}{2i-1} = \frac4\pi\sum_{i=1}^{N}\frac1{2i-1}$ 以下です。$i \ge 2$ では $y = 2/(2i-3)$ として $\log\frac{2i-1}{2i-3} = \log(1+y) \ge \frac{y}{1+y} = \frac{2}{2i-1}$ なので、 $$ \sum_{i=1}^{N}\frac1{2i-1} \le 1 + \frac12\sum_{i=2}^{N}\log\frac{2i-1}{2i-3} = 1 + \frac12\log(2N-1) $$ です。$2N - 1 = 2n+1$ なので $\lambda_n(x) \le 1 + \frac4\pi + \frac2\pi\log(2n+1)$ で、$x$ について最大をとれば主張を得ます。 $[a,b]$ の Chebyshev 点 $x_j = \frac{a+b}{2} + \frac{b-a}{2}\hat x_j$ では、$t(x) = (2x - a - b)/(b - a)$ について $\ell_j(x) = \hat\ell_j(t(x))$（$\hat\ell_j$ は $[-1,1]$ の点 $\hat x_j$ の Lagrange 基底）です。$x - x_k = \frac{b-a}{2}(t(x) - \hat x_k)$ の比をとれば因子 $\frac{b-a}{2}$ が約分されるからです。$t$ は $[a,b]$ を $[-1,1]$ の上へ写すので、Lebesgue 定数は等しくなります。