コンテンツにスキップ

多項式補間:誤差公式・Runge 現象・Chebyshev 点と Lebesgue 定数

対象読者: 学部 1〜2 年の微積分(平均値の定理・Rolle の定理・Taylor の定理)と、線形代数の行列式(Vandermonde 行列式を含む)を学んだ人。本テーマ第 1〜3 章の内容を前提にしてよく、補間の知識は仮定しません。

前提となる記事:平均値の定理とテイラーの定理:微分から関数の姿を復元する、行列式とその性質:符号付き体積としての det

生 Markdown

何点かでしか値がわからない関数を、多項式で置き換えて近似するにはどうすればよいか、という話です。点を通る多項式は簡単に作れますが、点をどう並べるかで結果が大きく変わり、等間隔に並べると点を増やすほどかえって悪くなる場合があります。この記事を読み終えると、補間多項式の作り方と誤差の見積もり方、そして点の置き方の良し悪しを数で測る方法がわかります。

  • 相異なる n+1n+1 点での値を与えると、それを通る nn 次以下の多項式(補間多項式)がただ 1 つ存在します。証明は Vandermonde 行列式が 00 でないことです。
  • 補間多項式は差分商を係数とする Newton 形で書け、点を 1 つ足すたびに項を 1 つ足すだけで更新できます。差分商は f(n)(ξ)/n!f^{(n)}(\xi)/n! に等しく、導関数の離散版です。
  • 補間誤差は f(x)−pn(x)=f(n+1)(ξ)(n+1)!∏j(x−xj)f(x) - p_n(x) = \dfrac{f^{(n+1)}(\xi)}{(n+1)!}\prod_j (x - x_j) です。証明は Rolle の定理の繰り返しです。値と導関数の両方を合わせる Hermite 補間も同じ形の誤差公式をもち、これが Gauss 求積の誤差評価の土台になります。
  • 補間の安定性は Lebesgue 定数 Λn\Lambda_n で測ります。∥f−pn∥≤(1+Λn)En(f)\|f - p_n\| \le (1 + \Lambda_n) E_n(f) で、Λn\Lambda_n はデータの誤差が何倍に増幅されるかを表す条件数でもあります。
  • 等間隔点では Λn≥2n/(4n2)\Lambda_n \ge 2^n/(4n^2) と指数的に増え、Runge の関数 1/(1+25x2)1/(1+25x^2) では補間が発散します。Chebyshev 点では Λn≤(2/π)log⁡(2n+1)+1+4/π\Lambda_n \le (2/\pi)\log(2n+1) + 1 + 4/\pi と対数的にしか増えず、解析関数なら誤差は幾何級数的に減ります。

関数 ff の値がいくつかの点でしかわからないとき、あるいは ff が複雑で何度も評価したくないとき、ff を多項式で置き換えて計算します。数値積分(多項式を積分する)、数値微分(多項式を微分する)、常微分方程式の多段法(過去の値を多項式で外挿する)、非線形方程式の割線法(2 点を通る 1 次式の零点)は、どれもこの置き換えの上に立っています。

自然な発想は「点を増やせば近似はよくなるだろう」というものです。ところが 1901 年に C. Runge が示したように、[−1,1][-1,1] の等間隔点で 1/(1+25x2)1/(1+25x^2) を補間すると、点を増やすほど区間の端の近くで誤差が爆発します。滑らかな(実解析的な)関数でも、点の配置を誤ると補間は収束しないのです。一方、区間の端に点を密に置く Chebyshev 点を使うと、同じ関数の補間は幾何級数的に収束します。

この違いを説明するのが本章の主題です。誤差公式の中の積 ∏j(x−xj)\prod_j(x - x_j) の大きさと、Lebesgue 定数という「補間という操作の条件数」の 2 つを調べると、等間隔点と Chebyshev 点の差が定量的にわかります。

以下、Pn\mathcal{P}_n で実係数の nn 次以下の多項式全体を表します。

補題 2.1(Vandermonde 行列式)✓

n≥1n \ge 1 とし、x0,…,xnx_0, \ldots, x_n を実数とする。(i,j)(i, j) 成分が xijx_i^j(0≤i,j≤n0 \le i, j \le n)の行列 VV について

det⁡V=∏0≤i<j≤n(xj−xi)\det V = \prod_{0 \le i \lt j \le n} (x_j - x_i)

が成り立つ。とくに x0,…,xnx_0, \ldots, x_n が相異なれば det⁡V≠0\det V \ne 0 である。

証明(補題 2.1)

nn についての帰納法で示します。n=1n = 1 では det⁡V=x1−x0\det V = x_1 - x_0 で正しいです。n≥2n \ge 2 とし、nn 個の点についての主張を仮定します。

VV の第 jj 列から第 j−1j-1 列の x0x_0 倍を引く操作を、j=n,n−1,…,1j = n, n-1, \ldots, 1 の順に行います(後ろから行うので、引く列はつねに元のままです)。この操作で行列式は変わりません(命題 4.3[行列式とその性質] の 3)。操作後の (i,j)(i, j) 成分は、j≥1j \ge 1 で xij−x0xij−1=xij−1(xi−x0)x_i^j - x_0 x_i^{j-1} = x_i^{j-1}(x_i - x_0) です。とくに第 00 行は (1,0,…,0)(1, 0, \ldots, 0) になります。第 00 行について余因子展開(定理 5.2[行列式とその性質])すると、det⁡V\det V は第 ii 行(i=1,…,ni = 1, \ldots, n)が (xi−x0)(1,xi,…,xin−1)(x_i - x_0)(1, x_i, \ldots, x_i^{n-1}) の nn 次の行列式に等しくなります。各行から xi−x0x_i - x_0 をくくり出すと(命題 4.3[行列式とその性質] の 2 を、5 により行に適用)

det⁡V=∏i=1n(xi−x0)⋅det⁡(xij)1≤i≤n, 0≤j≤n−1\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}

です。右辺の行列式は x1,…,xnx_1, \ldots, x_n の Vandermonde 行列式なので、帰納法の仮定から ∏1≤i<j≤n(xj−xi)\prod_{1 \le i \lt j \le n}(x_j - x_i) に等しく、主張を得ます。

定理 2.2(補間多項式の存在と一意性)✓

x0,…,xnx_0, \ldots, x_n を相異なる実数、y0,…,yny_0, \ldots, y_n を実数とする。p(xi)=yip(x_i) = y_i(i=0,…,ni = 0, \ldots, n)を満たす p∈Pnp \in \mathcal{P}_n がただ 1 つ存在し、

p(x)=∑i=0nyi ℓi(x),ℓi(x)=∏j≠ix−xjxi−xjp(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 形)。

証明(定理 2.2)

n=0n = 0 なら定数 y0y_0 がただ 1 つの解です。n≥1n \ge 1 とします。p(x)=∑j=0najxjp(x) = \sum_{j=0}^n a_j x^j とおくと、条件 p(xi)=yip(x_i) = y_i は係数ベクトル a=(a0,…,an)Ta = (a_0, \ldots, a_n)^{\mathsf{T}} についての連立 1 次方程式 Va=yVa = y です(VV は 補題 2.1 の行列)。det⁡V≠0\det V \ne 0 なので、正則性の判定(定理 7.1[行列式とその性質] の 1 ⇒ 5)からただ 1 つの解があります。ℓi\ell_i は nn 次の多項式で ℓi(xk)=1\ell_i(x_k) = 1(k=ik = i)、00(k≠ik \ne i)を満たすので、∑iyiℓi\sum_i y_i\ell_i は条件を満たす Pn\mathcal{P}_n の元であり、一意性からこれが解です。

系 2.3(零点の個数)✓

q∈Pnq \in \mathcal{P}_n が相異なる n+1n + 1 個の点で 00 になるならば、qq は零多項式である。

証明(系 2.3)

qq と零多項式は、どちらもその n+1n+1 点で値 00 をとる Pn\mathcal{P}_n の元です。定理 2.2 の一意性から両者は等しくなります。

以下、ff の値 yi=f(xi)y_i = f(x_i) を補間する Pn\mathcal{P}_n の元を pnp_n(点を明示したいときは px0,…,xnp_{x_0, \ldots, x_n})と書き、ω(x)=∏j=0n(x−xj)\omega(x) = \prod_{j=0}^{n}(x - x_j) とおきます。Lagrange 形は理論には便利ですが、点を 1 つ足すとすべての ℓi\ell_i が変わります。次の節の Newton 形はこの欠点を解消します。

定義 3.1(差分商)↗

相異なる点 x0,…,xkx_0, \ldots, x_k と関数 ff について、ff をこの k+1k+1 点で補間する Pk\mathcal{P}_k の元の xkx^k の係数を f[x0,…,xk]f[x_0, \ldots, x_k] と書き、ff の kk 階の差分商という。f[x0]=f(x0)f[x_0] = f(x_0) である。

命題 3.2(差分商の対称性と陽な表示)✓

f[x0,…,xk]f[x_0, \ldots, x_k] は点の並べ方によらず、

f[x0,…,xk]=∑i=0kf(xi)∏j≠i(xi−xj)f[x_0, \ldots, x_k] = \sum_{i=0}^{k}\frac{f(x_i)}{\prod_{j \ne i}(x_i - x_j)}

である。

証明(命題 3.2)

補間多項式は点の集合だけで決まる(定理 2.2)ので、並べ方によりません。Lagrange 形の ℓi\ell_i の xkx^k の係数は 1/∏j≠i(xi−xj)1/\prod_{j \ne i}(x_i - x_j) なので、∑if(xi)ℓi\sum_i f(x_i)\ell_i の xkx^k の係数は右辺です。

定理 3.3(Newton 形と差分商の漸化式)✓

x0,…,xnx_0, \ldots, x_n を相異なる点とする。

  1. ff を補間する pn∈Pnp_n \in \mathcal{P}_n は pn(x)=∑k=0nf[x0,…,xk]∏j=0k−1(x−xj)p_n(x) = \sum_{k=0}^{n} f[x_0, \ldots, x_k]\prod_{j=0}^{k-1}(x - x_j) と書ける(Newton 形。k=0k = 0 の積は 11)。
  2. k≥1k \ge 1 について f[x0,…,xk]=f[x1,…,xk]−f[x0,…,xk−1]xk−x0f[x_0, \ldots, x_k] = \frac{f[x_1, \ldots, x_k] - f[x_0, \ldots, x_{k-1}]}{x_k - x_0} が成り立つ。
証明(定理 3.3)

1. pk=px0,…,xkp_k = p_{x_0, \ldots, x_k} とおき、pk=pk−1+f[x0,…,xk]∏j<k(x−xj)p_k = p_{k-1} + f[x_0, \ldots, x_k]\prod_{j \lt k}(x - x_j)(k≥1k \ge 1)を示せば、kk について足し合わせて主張を得ます。r=pk−pk−1r = p_k - p_{k-1} は Pk\mathcal{P}_k の元で x0,…,xk−1x_0, \ldots, x_{k-1} で 00 になり、pk−1p_{k-1} の次数は k−1k-1 以下なので xkx^k の係数は f[x0,…,xk]f[x_0, \ldots, x_k] です。s=r−f[x0,…,xk]∏j<k(x−xj)s = r - f[x_0, \ldots, x_k]\prod_{j \lt k}(x - x_j) は xkx^k の項が打ち消し合うので Pk−1\mathcal{P}_{k-1} の元で、kk 個の点 x0,…,xk−1x_0, \ldots, x_{k-1} で 00 になります。系 2.3 から s=0s = 0 です。

2. q0=px0,…,xk−1q_0 = p_{x_0, \ldots, x_{k-1}}、q1=px1,…,xkq_1 = p_{x_1, \ldots, x_k}(ともに Pk−1\mathcal{P}_{k-1} の元)とし、

P(x)=(x−x0) q1(x)−(x−xk) q0(x)xk−x0P(x) = \frac{(x - x_0)\,q_1(x) - (x - x_k)\,q_0(x)}{x_k - x_0}

とおきます。P∈PkP \in \mathcal{P}_k で、P(x0)=q0(x0)=f(x0)P(x_0) = q_0(x_0) = f(x_0)、P(xk)=q1(xk)=f(xk)P(x_k) = q_1(x_k) = f(x_k)、0<i<k0 \lt i \lt k では P(xi)=((xi−x0)−(xi−xk))f(xi)/(xk−x0)=f(xi)P(x_i) = \bigl((x_i - x_0) - (x_i - x_k)\bigr)f(x_i)/(x_k - x_0) = f(x_i) です。一意性(定理 2.2)から P=pkP = p_k で、両辺の xkx^k の係数を比べると漸化式を得ます。

漸化式 2 から、差分商は表の形で次々に計算できます(f[xi,xi+1]f[x_i, x_{i+1}] の列、f[xi,xi+1,xi+2]f[x_i, x_{i+1}, x_{i+2}] の列、…)。k=1,2k = 1, 2 の場合は、非線形方程式の反復解法で割線法の解析に使った差分商と一致します。そこで示した 2 階の場合(補題 6.1[非線形方程式の反復解法])は、定理 3.5 の n=2n = 2 の場合です。

補題 3.4(Rolle の定理の繰り返し)✓

II を開区間、k≥1k \ge 1 とし、g:I→Rg : I \to \mathbb{R} は kk 回微分可能で、相異なる k+1k + 1 個の点 t0<t1<⋯<tkt_0 \lt t_1 \lt \cdots \lt t_k(II の点)で 00 になるとする。このとき g(k)(ξ)=0g^{(k)}(\xi) = 0 を満たす ξ∈(t0,tk)\xi \in (t_0, t_k) が存在する。

証明(補題 3.4)

kk についての帰納法です。gg は微分可能なので各 [ti,ti+1][t_{i}, t_{i+1}] で連続かつ内部で微分可能で、両端で 00 なので、Rolle の定理(定理 3.1[平均値の定理とテイラーの定理])から g′(si)=0g'(s_i) = 0 となる si∈(ti,ti+1)s_i \in (t_i, t_{i+1}) があります(i=0,…,k−1i = 0, \ldots, k-1)。s0<⋯<sk−1s_0 \lt \cdots \lt s_{k-1} は相異なり (t0,tk)(t_0, t_k) に入ります。k=1k = 1 ならこれで終わりです。k≥2k \ge 2 なら g′g' は k−1k-1 回微分可能で kk 個の点で 00 なので、帰納法の仮定から (g′)(k−1)=g(k)(g')^{(k-1)} = g^{(k)} は (s0,sk−1)⊂(t0,tk)(s_0, s_{k-1}) \subset (t_0, t_k) の点で 00 になります。

定理 3.5(差分商の平均値定理)✓

n≥1n \ge 1、II を開区間、f:I→Rf : I \to \mathbb{R} を nn 回微分可能とし、x0,…,xn∈Ix_0, \ldots, x_n \in I を相異なる点とする。このとき min⁡ixi<ξ<max⁡ixi\min_i x_i \lt \xi \lt \max_i x_i を満たす ξ\xi があって

f[x0,…,xn]=f(n)(ξ)n!f[x_0, \ldots, x_n] = \frac{f^{(n)}(\xi)}{n!}

が成り立つ。

証明(定理 3.5)

g=f−png = f - p_n は II 上 nn 回微分可能で、n+1n+1 個の点 x0,…,xnx_0, \ldots, x_n で 00 です。補題 3.4 から g(n)(ξ)=0g^{(n)}(\xi) = 0 となる ξ\xi が最小の点と最大の点の間にあります。pnp_n の xnx^n の係数は f[x0,…,xn]f[x_0, \ldots, x_n] なので pn(n)=n! f[x0,…,xn]p_n^{(n)} = n!\,f[x_0, \ldots, x_n] で、g(n)(ξ)=0g^{(n)}(\xi) = 0 が主張です。

定理 4.1(補間誤差の公式)✓

II を開区間、f:I→Rf : I \to \mathbb{R} を n+1n+1 回微分可能とし、x0,…,xn∈Ix_0, \ldots, x_n \in I を相異なる点、pn∈Pnp_n \in \mathcal{P}_n をそこでの ff の補間多項式とする。任意の x∈Ix \in I について、x,x0,…,xnx, x_0, \ldots, x_n の最小値と最大値の間の ξ\xi があって

f(x)−pn(x)=f(n+1)(ξ)(n+1)! ω(x),ω(x)=∏j=0n(x−xj)f(x) - p_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\,\omega(x), \qquad \omega(x) = \prod_{j=0}^{n}(x - x_j)

が成り立つ。

証明の方針は、誤差を ω(x)\omega(x) の定数倍で打ち消す補助関数 gg を作り、gg の零点が多いことから Rolle の定理(補題 3.4)で g(n+1)g^{(n+1)} の零点を取り出す、というものです。なお「モニック」は最高次の係数が 11 の多項式のことで、ω\omega はその例です。

証明(定理 4.1)

xx がどれかの xjx_j に等しければ両辺は 00 です。そうでないとき ω(x)≠0\omega(x) \ne 0 なので、c=(f(x)−pn(x))/ω(x)c = (f(x) - p_n(x))/\omega(x) とおき、g(t)=f(t)−pn(t)−c ω(t)g(t) = f(t) - p_n(t) - c\,\omega(t) を考えます。gg は n+1n+1 回微分可能で、n+2n+2 個の相異なる点 x0,…,xn,xx_0, \ldots, x_n, x で 00 です。補題 3.4(k=n+1k = n+1)から、それらの最小値と最大値の間の ξ\xi で g(n+1)(ξ)=0g^{(n+1)}(\xi) = 0 です。pn(n+1)=0p_n^{(n+1)} = 0、ω\omega はモニックな n+1n+1 次式なので ω(n+1)=(n+1)!\omega^{(n+1)} = (n+1)! で、0=f(n+1)(ξ)−c(n+1)!0 = f^{(n+1)}(\xi) - c(n+1)! です。cc の定義に戻せば主張を得ます。

誤差は 2 つの因子の積です。f(n+1)/(n+1)!f^{(n+1)}/(n+1)! は関数の性質で、こちらは選べません。ω(x)\omega(x) は点の配置だけで決まり、こちらは選べます。区間の上で max⁡∣ω∣\max|\omega| を最小にする点の選び方が、後の節の Chebyshev 点です(系 6.3)。

定理 4.2(Hermite 補間の存在・一意性と誤差)✓

x0,…,xnx_0, \ldots, x_n を相異なる実数とする。

  1. 任意の実数 yi,ziy_i, z_i(i=0,…,ni = 0, \ldots, n)について、H(xi)=yiH(x_i) = y_i、H′(xi)=ziH'(x_i) = z_i を満たす H∈P2n+1H \in \mathcal{P}_{2n+1} がただ 1 つ存在する。
  2. II を x0,…,xnx_0, \ldots, x_n を含む開区間、f:I→Rf : I \to \mathbb{R} を 2n+22n+2 回微分可能とし、HH を yi=f(xi)y_i = f(x_i)、zi=f′(xi)z_i = f'(x_i) についての 1 の多項式とする。任意の x∈Ix \in I について、x,x0,…,xnx, x_0, \ldots, x_n の最小値と最大値の間の ξ\xi があって f(x)−H(x)=f(2n+2)(ξ)(2n+2)! ω(x)2f(x) - H(x) = \frac{f^{(2n+2)}(\xi)}{(2n+2)!}\,\omega(x)^2 が成り立つ。

証明の方針は、1 では未知数と条件の数がどちらも 2n+22n+2 なので、すべて 00 のデータから H=0H = 0 しか出ないことを示せば十分だ、というものです。2 は補間誤差の公式と同じ手順で、各 xix_i が値と傾きの 2 回分の零点として数えられるため、零点が 2n+22n+2 個になって ω(x)2\omega(x)^2 が現れます。

証明(定理 4.2)

1. H(x)=∑k=02n+1akxkH(x) = \sum_{k=0}^{2n+1} a_k x^k とおくと、2n+22n+2 個の条件は係数 a0,…,a2n+1a_0, \ldots, a_{2n+1} についての 2n+22n+2 元の連立 1 次方程式です。正則性の判定(定理 7.1[行列式とその性質] の 4 ⇒ 5)から、斉次方程式(yi=zi=0y_i = z_i = 0)の解が H=0H = 0 だけであることを示せば十分です。H∈P2n+1H \in \mathcal{P}_{2n+1} が H(xi)=H′(xi)=0H(x_i) = H'(x_i) = 0 を満たすとします。点を小さい順に並べると、Rolle の定理(定理 3.1[平均値の定理とテイラーの定理])から隣り合う 2 点の間に H′H' の零点が 1 つずつ、計 nn 個あり、これらは xix_i と異なります。H′(xi)=0H'(x_i) = 0 と合わせて、H′∈P2nH' \in \mathcal{P}_{2n} は相異なる 2n+12n+1 個の零点をもつので、系 2.3 から H′=0H' = 0 です。よって HH は定数で、H(x0)=0H(x_0) = 0 から H=0H = 0 です。

2. xx がどれかの xjx_j なら両辺は 00 です。そうでないとき c=(f(x)−H(x))/ω(x)2c = (f(x) - H(x))/\omega(x)^2 とおき、g(t)=f(t)−H(t)−c ω(t)2g(t) = f(t) - H(t) - c\,\omega(t)^2 を考えます。gg は n+2n + 2 個の相異なる点 x0,…,xn,xx_0, \ldots, x_n, x で 00 なので、Rolle の定理から、これらを小さい順に並べたときの隣り合う 2 点の間に g′g' の零点が 1 つずつ、計 n+1n+1 個あります。さらに ω(xi)=0\omega(x_i) = 0 から (ω2)′(xi)=2ω(xi)ω′(xi)=0(\omega^2)'(x_i) = 2\omega(x_i)\omega'(x_i) = 0 なので g′(xi)=f′(xi)−H′(xi)=0g'(x_i) = f'(x_i) - H'(x_i) = 0 です。xix_i は並べた点そのものなので、点の間にある先の n+1n+1 個とは異なります。よって g′g' は相異なる 2n+22n+2 個の零点をもち、補題 3.4(g′g' に k=2n+1k = 2n+1 で適用)から g(2n+2)(ξ)=0g^{(2n+2)}(\xi) = 0 となる ξ\xi が、x,x0,…,xnx, x_0, \ldots, x_n の最小値と最大値の間にあります。H(2n+2)=0H^{(2n+2)} = 0、ω2\omega^2 はモニックな 2n+22n+2 次式なので (ω2)(2n+2)=(2n+2)!(\omega^2)^{(2n+2)} = (2n+2)! で、f(2n+2)(ξ)=c (2n+2)!f^{(2n+2)}(\xi) = c\,(2n+2)! です。

Hermite 補間は、n+1n+1 点の値と傾きという 2n+22n+2 個の情報から 2n+12n+1 次式を決めます。誤差公式の ω(x)2\omega(x)^2 は区間の上で符号を変えないので、積分の平均値の定理と組み合わせると「f(2n+2)(ξ)/(2n+2)!f^{(2n+2)}(\xi)/(2n+2)! × ∫ω2\int\omega^2」の形の誤差公式が得られます。これが Newton–Cotes 則と Gauss 求積の章で Gauss 求積の誤差を導く道筋です。

例 4.3(Runge の関数を等間隔点と Chebyshev 点で補間する)

f(x)=1/(1+25x2)f(x) = 1/(1+25x^2) を [−1,1][-1,1] の n+1n+1 個の等間隔点と、Tn+1T_{n+1} の零点 cos⁡((2k+1)π/(2n+2))\cos\bigl((2k+1)\pi/(2n+2)\bigr)(後の節で定義する Chebyshev 点)で補間し、4001 点の格子上の最大誤差を比べます。補間多項式の係数は求めず、SciPy の重心形の補間(注意 6.9)で値を計算します。

# 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}")
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≥5n \ge 5 の補間は 5 次式を正確に再現するはずで(定理 2.2 の一意性)、最後の列は誤差が 10−1210^{-12} 未満であることを示しています。(b) Chebyshev 点の誤差は、重心形の計算と、NumPy の chebinterpolate(同じ点での補間を Chebyshev 級数の係数として求める別の算法)とで表示の桁まで一致します。(c) 等間隔点の誤差は nn を 2 倍にするごとに大きくなり(Runge 現象)、Chebyshev 点の誤差は nn を 10 増やすごとにほぼ一定の比で減ります。この比が ρ−10\rho^{-10}(ρ=(1+26)/5\rho = (1+\sqrt{26})/5)であることは、後の例(例 6.8)で確かめます。

n=5n = 5 では等間隔点のほうがわずかに良いことにも注意してください。Chebyshev 点が有利なのは点数を増やしたときの振る舞いで、どの nn でも勝つわけではありません。等間隔点で何が起きているのかを、次の節で補間という操作そのものの性質として捉えます。

この節では、補間という操作が入力のずれをどれだけ増幅するかを 1 つの数で測ります。[a,b][a, b] 上の連続関数の空間 C[a,b]C[a,b] に一様ノルム ∥g∥=max⁡[a,b]∣g∣\|g\| = \max_{[a,b]}|g|(関数の絶対値の最大値)を入れます。[a,b][a,b] の相異なる点 x0,…,xnx_0, \ldots, x_n を固定し、f∈C[a,b]f \in C[a,b] にその補間多項式を対応させる写像を Inf=∑if(xi)ℓi\mathcal{I}_n f = \sum_i f(x_i)\ell_i と書きます。In\mathcal{I}_n は線形で、定理 2.2 の一意性から q∈Pnq \in \mathcal{P}_n について Inq=q\mathcal{I}_n q = q です。

定義 5.1(Lebesgue 関数と Lebesgue 定数)↗

λn(x)=∑i=0n∣ℓi(x)∣\lambda_n(x) = \sum_{i=0}^{n}|\ell_i(x)| を点 x0,…,xnx_0, \ldots, x_n の Lebesgue 関数、Λn=max⁡x∈[a,b]λn(x)\Lambda_n = \max_{x \in [a,b]}\lambda_n(x) を Lebesgue 定数という。また f∈C[a,b]f \in C[a,b] について En(f)=inf⁡q∈Pn∥f−q∥E_n(f) = \inf_{q \in \mathcal{P}_n}\|f - q\| を、nn 次以下の多項式による最良一様近似の誤差という。

λn\lambda_n は連続なので、最大値は存在します(定理 2.3[平均値の定理とテイラーの定理])。∑iℓi=1\sum_i \ell_i = 1(定数 11 の補間)なので λn≥1\lambda_n \ge 1、したがって Λn≥1\Lambda_n \ge 1 です。

次の定理の 1 の「作用素ノルム」は、∥f∥=1\|f\| = 1 の関数 ff に補間を施したとき、結果の大きさが最大でどこまで大きくなるかを表す数のことです。1 はそれが Λn\Lambda_n に等しいと述べていて、証明の後半は、その最大を実現する折れ線の関数を実際に作る部分です。

定理 5.2(Lebesgue 定数による誤差の評価)✓

[a,b][a,b] の相異なる点 x0,…,xnx_0, \ldots, x_n と f∈C[a,b]f \in C[a,b] について次が成り立つ。

  1. ∥Inf∥≤Λn∥f∥\|\mathcal{I}_n f\| \le \Lambda_n\|f\| であり、等号を満たす f≠0f \ne 0 が存在する。すなわち Λn\Lambda_n は In\mathcal{I}_n の作用素ノルムである。
  2. ∥f−Inf∥≤(1+Λn) En(f)\|f - \mathcal{I}_n f\| \le (1 + \Lambda_n)\,E_n(f)。
  3. En(f)→0E_n(f) \to 0(n→∞n \to \infty)。
証明(定理 5.2)

1. ∣Inf(x)∣≤∑i∣f(xi)∣∣ℓi(x)∣≤∥f∥λn(x)≤∥f∥Λn|\mathcal{I}_n f(x)| \le \sum_i |f(x_i)||\ell_i(x)| \le \|f\|\lambda_n(x) \le \|f\|\Lambda_n です。等号のために、λn(x^)=Λn\lambda_n(\hat x) = \Lambda_n となる x^\hat x を取り、sis_i を ℓi(x^)\ell_i(\hat x) の符号(00 なら 11)とします。点を小さい順に並べ、ff を各点で sis_i をとる折れ線(両端の点の外では一定)とすると、ff は連続で ∥f∥=1\|f\| = 1 であり、Inf(x^)=∑isiℓi(x^)=λn(x^)=Λn\mathcal{I}_n f(\hat x) = \sum_i s_i\ell_i(\hat x) = \lambda_n(\hat x) = \Lambda_n です。

2. 任意の q∈Pnq \in \mathcal{P}_n について Inq=q\mathcal{I}_n q = q なので、f−Inf=(f−q)−In(f−q)f - \mathcal{I}_n f = (f - q) - \mathcal{I}_n(f - q) です。1 から ∥f−Inf∥≤(1+Λn)∥f−q∥\|f - \mathcal{I}_n f\| \le (1 + \Lambda_n)\|f - q\| で、qq について下限をとれば主張を得ます。

3. ε>0\varepsilon \gt 0 について、Weierstrass の多項式近似定理(系 4.4[Arzelà–Ascoli と Stone–Weierstrass])から ∥f−P∥<ε\|f - P\| \lt \varepsilon となる多項式 PP があります。PP の次数を NN とすると、n≥Nn \ge N で P∈PnP \in \mathcal{P}_n なので En(f)≤εE_n(f) \le \varepsilon です。

2 と 3 を合わせても、∥f−Inf∥→0\|f - \mathcal{I}_n f\| \to 0 とは言えません。En(f)E_n(f) が 00 に近づく速さより Λn\Lambda_n が増える速さのほうが速ければ、補間は発散しえます。Runge の現象はまさにこの状況です。

例 5.3(Lebesgue 定数は補間の条件数)

データ f(xi)f(x_i) に大きさ ε\varepsilon 以下の誤差(測定誤差や丸め誤差)δi\delta_i が乗ると、補間多項式は ∑iδiℓi\sum_i\delta_i\ell_i だけずれます。その一様ノルムは εΛn\varepsilon\Lambda_n 以下で、定理 5.2 の 1 から、最悪の場合(δi=εsi\delta_i = \varepsilon s_i)にはちょうど εΛn\varepsilon\Lambda_n になります。つまり Λn\Lambda_n は「データから補間多項式への写像」の(一様ノルムで測った絶対)条件数です。倍精度の丸め誤差 u≈1.1×10−16u \approx 1.1 \times 10^{-16} を想定すると、Λn≈1010\Lambda_n \approx 10^{10} の点の配置では補間の値は最悪で 6 桁程度しか信用できません。後の表(例 6.6)では等間隔点の Λ30\Lambda_{30} が約 6.6×1066.6 \times 10^6 で、nn を 1 増やすごとにおよそ 2 倍になっていくので、n=40n = 40 では約 5×1095 \times 10^9 とこの水準に迫ります。

命題 5.4(等間隔点の Lebesgue 定数の下界)✓

n≥1n \ge 1、h=(b−a)/nh = (b-a)/n、xj=a+jhx_j = a + jh(j=0,…,nj = 0, \ldots, n)とする。このとき

Λn≥λn(a+h2)≥2n4n2\Lambda_n \ge \lambda_n\Bigl(a + \frac h2\Bigr) \ge \frac{2^n}{4n^2}

である。

証明の方針は、端の最初の小区間の中点 a+h/2a + h/2 だけに注目し、そこでの各 ∣ℓj∣|\ell_j| を下から評価して足し合わせる、というものです。足し合わせには二項定理を使うので、式変形が続きますが、最後は (2nn)\binom{2n}{n} の下からの評価に帰着します。

証明(命題 5.4)

x=a+thx = a + th とおくと ℓj(x)=∏k≠j(t−k)/(j−k)\ell_j(x) = \prod_{k \ne j}(t - k)/(j - k) です。t=1/2t = 1/2 では

∏k=0n∣12−k∣=12∏k=1n2k−12=(2n−1)!!2n+1,∏k≠j∣j−k∣=j! (n−j)!\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⋅3⋯(2n−1)(2n-1)!! = 1 \cdot 3 \cdots (2n-1))、

∣ℓj(a+h2)∣=(2n−1)!!2n+1 ∣1/2−j∣  j! (n−j)!≥(2n−1)!!2n+1 n  j! (n−j)!\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∣≤n|1/2 - j| \le n を使いました)。二項定理から ∑j=0n1/(j!(n−j)!)=2n/n!\sum_{j=0}^{n} 1/(j!(n-j)!) = 2^n/n! なので

λn(a+h2)≥(2n−1)!!2n⋅n!=(2n)!2n+1 n (n!)2=12n+1n(2nn)\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)!/(2nn!)(2n-1)!! = (2n)!/(2^n n!))。最後に (2nn)≥4n/(2n)\binom{2n}{n} \ge 4^n/(2n) を示します。(2nk+1)/(2nk)=(2n−k)/(k+1)\binom{2n}{k+1}/\binom{2n}{k} = (2n-k)/(k+1) は k≤n−1k \le n-1 で 11 以上、k≥nk \ge n で 11 未満なので、(2nn)\binom{2n}{n} は (2nk)\binom{2n}{k} の最大です。二項定理から 4n=∑k=02n(2nk)=2+∑k=12n−1(2nk)≤2+(2n−1)(2nn)≤2n(2nn)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}((2nn)≥2\binom{2n}{n} \ge 2 を使いました)です。合わせて λn(a+h/2)≥4n/(2n+1n⋅2n)=2n/(4n2)\lambda_n(a + h/2) \ge 4^n/(2^{n+1}n \cdot 2n) = 2^n/(4n^2) です。

下界は n=10n = 10 で 2.562.56、n=30n = 30 で約 3.0×1053.0 \times 10^5 と控えめですが、指数的に増えることが要点です。実際の値はもっと大きく(後の表で n=30n = 30 で約 6.6×1066.6 \times 10^6)、およそ 2n+1/(e nlog⁡n)2^{n+1}/(e\,n\log n) の速さで増えることが知られています(Trefethen『Approximation Theory and Approximation Practice』第 15 章)。

Chebyshev 多項式 TkT_k は、漸化式 T0=1T_0 = 1、T1=xT_1 = x、Tk+1=2xTk−Tk−1T_{k+1} = 2xT_k - T_{k-1} で定めます(定義 2.1[誤差保証つきの求積と関数近似])。

この多項式の要点は、x=cos⁡θx = \cos\theta とおくと Tk(x)=cos⁡kθT_k(x) = \cos k\theta になることです。角度 θ\theta を等間隔に動かすと cos⁡kθ\cos k\theta は等しい振れ幅で振動し、[−1,1][-1,1] 上で ∣Tk∣≤1|T_k| \le 1 が保たれます。零点も角度が等間隔の点に対応するので、Chebyshev 点は、半円周上に等間隔に並べた点を横軸へ落とした点になります。下の図は n=5n = 5(6 点)の場合です。

-101
半円周上で角度が等間隔の 6 点(15°, 45°, ..., 165°)を横軸へ真下に落とした点が、Chebyshev 点 cos((2j+1)π/12)。区間の端ほど点が密になる。

補題 6.1(Chebyshev 多項式の三角関数による表示)✓

k≥0k \ge 0 について次が成り立つ。

  1. すべての実数 θ\theta について Tk(cos⁡θ)=cos⁡kθT_k(\cos\theta) = \cos k\theta。とくに [−1,1][-1,1] 上で ∣Tk∣≤1|T_k| \le 1。
  2. k≥1k \ge 1 なら TkT_k は kk 次で、最高次の係数は 2k−12^{k-1} である。Tk(−x)=(−1)kTk(x)T_k(-x) = (-1)^kT_k(x)。
  3. k≥1k \ge 1 なら TkT_k の零点は cos⁡((2j+1)π/(2k))\cos\bigl((2j+1)\pi/(2k)\bigr)(j=0,…,k−1j = 0, \ldots, k-1)の kk 個で、すべて (−1,1)(-1, 1) の相異なる点である。また ηj=cos⁡(jπ/k)\eta_j = \cos(j\pi/k)(j=0,…,kj = 0, \ldots, k)で Tk(ηj)=(−1)jT_k(\eta_j) = (-1)^j である。
証明(補題 6.1)

1. k=0,1k = 0, 1 では定義どおりです。k−1,kk - 1, k で成り立てば、加法定理 cos⁡(k+1)θ+cos⁡(k−1)θ=2cos⁡θcos⁡kθ\cos(k+1)\theta + \cos(k-1)\theta = 2\cos\theta\cos k\theta から

Tk+1(cos⁡θ)=2cos⁡θ Tk(cos⁡θ)−Tk−1(cos⁡θ)=2cos⁡θcos⁡kθ−cos⁡(k−1)θ=cos⁡(k+1)θ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∈[−1,1]x \in [-1,1] は x=cos⁡θx = \cos\theta と書けるので ∣Tk(x)∣≤1|T_k(x)| \le 1 です。

2. 漸化式から帰納法で、Tk+1T_{k+1} の最高次の項は 2x⋅2k−1xk=2kxk+12x \cdot 2^{k-1}x^k = 2^k x^{k+1} です(k=0k = 0 から 11 へは T1=xT_1 = x)。偶奇も漸化式から帰納法で従います。

3. θj=(2j+1)π/(2k)∈(0,π)\theta_j = (2j+1)\pi/(2k) \in (0, \pi) は相異なり、cos⁡\cos は [0,π][0, \pi] で狭義減少なので cos⁡θj\cos\theta_j は (−1,1)(-1,1) の相異なる kk 点です。1 から Tk(cos⁡θj)=cos⁡((2j+1)π/2)=0T_k(\cos\theta_j) = \cos\bigl((2j+1)\pi/2\bigr) = 0 です。ほかにも零点があれば、Tk∈PkT_k \in \mathcal{P}_k が k+1k+1 個の相異なる零点をもつので 系 2.3 から Tk=0T_k = 0 となり、最高次の係数が 2k−1≠02^{k-1} \ne 0 であることに矛盾します。よってこれがすべての零点です。Tk(cos⁡(jπ/k))=cos⁡jπ=(−1)jT_k(\cos(j\pi/k)) = \cos j\pi = (-1)^j も 1 から従います。

定理 6.2(Chebyshev 多項式の最小最大性)✓

n≥0n \ge 0 とする。

  1. [−1,1][-1,1] 上で、n+1n+1 次のモニック多項式 qq は max⁡[−1,1]∣q∣≥2−n\max_{[-1,1]}|q| \ge 2^{-n} を満たし、q=2−nTn+1q = 2^{-n}T_{n+1} で等号が成り立つ。
  2. a<ba \lt b とする。[a,b][a,b] 上で、n+1n+1 次のモニック多項式 qq は max⁡[a,b]∣q∣≥2((b−a)/4)n+1\max_{[a,b]}|q| \ge 2\bigl((b-a)/4\bigr)^{n+1} を満たす。
証明(定理 6.2)

1. T~=2−nTn+1\tilde T = 2^{-n}T_{n+1} は 補題 6.1 の 2 からモニックで、1 から max⁡∣T~∣≤2−n\max|\tilde T| \le 2^{-n}、3 から ηj=cos⁡(jπ/(n+1))\eta_j = \cos(j\pi/(n+1))(j=0,…,n+1j = 0, \ldots, n+1)で T~(ηj)=(−1)j2−n\tilde T(\eta_j) = (-1)^j2^{-n} です。よって max⁡∣T~∣=2−n\max|\tilde T| = 2^{-n} です。

max⁡∣q∣<2−n\max|q| \lt 2^{-n} を満たすモニックな qq があったと仮定します。r=T~−qr = \tilde T - q は最高次の項が打ち消し合うので Pn\mathcal{P}_n の元です。∣q(ηj)∣<2−n|q(\eta_j)| \lt 2^{-n} なので r(ηj)=(−1)j2−n−q(ηj)r(\eta_j) = (-1)^j2^{-n} - q(\eta_j) の符号は (−1)j(-1)^j です。η0>η1>⋯>ηn+1\eta_0 \gt \eta_1 \gt \cdots \gt \eta_{n+1} の隣り合う 2 点で rr の符号は逆なので、中間値の定理(定理 5.7[極限と連続性])から各区間 (ηj+1,ηj)(\eta_{j+1}, \eta_j) に rr の零点があり、相異なる零点が n+1n+1 個得られます。系 2.3 から r=0r = 0 ですが、r(η0)>0r(\eta_0) \gt 0 に矛盾します。

2. x=a+b2+b−a2tx = \frac{a+b}{2} + \frac{b-a}{2}t とおくと、Q(t)=(2/(b−a))n+1q(a+b2+b−a2t)Q(t) = (2/(b-a))^{n+1}q\bigl(\frac{a+b}{2} + \frac{b-a}{2}t\bigr) は tt の n+1n+1 次のモニック多項式です。1 から max⁡t∈[−1,1]∣Q(t)∣≥2−n\max_{t \in [-1,1]}|Q(t)| \ge 2^{-n} なので、max⁡[a,b]∣q∣≥((b−a)/2)n+12−n=2((b−a)/4)n+1\max_{[a,b]}|q| \ge ((b-a)/2)^{n+1}2^{-n} = 2((b-a)/4)^{n+1} です。

この定理を補間誤差の公式に組み合わせると、max⁡∣ω∣\max|\omega| を最小にする点が決まります。

系 6.3(Chebyshev 点での補間誤差)✓

a<ba \lt b、n≥0n \ge 0 とし、Chebyshev 点

xj=a+b2+b−a2cos⁡(2j+1)π2n+2(j=0,…,n)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]∣ω∣=2((b−a)/4)n+1\max_{[a,b]}|\omega| = 2((b-a)/4)^{n+1} で、これは相異なる n+1n+1 点のすべての選び方の中で最小である。ff が [a,b][a,b] を含む開区間で n+1n+1 回微分可能で ∣f(n+1)∣≤Mn+1|f^{(n+1)}| \le M_{n+1}([a,b][a,b] 上)ならば

max⁡[a,b]∣f−pn∣≤2Mn+1(n+1)!(b−a4)n+1\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][a,b] = [-1,1] では右辺は Mn+1/(2n(n+1)!)M_{n+1}/(2^n(n+1)!) である。

証明(系 6.3)

t(x)=(2x−a−b)/(b−a)t(x) = (2x - a - b)/(b-a) とおきます。補題 6.1 の 2・3 から、xx の多項式 Tn+1(t(x))T_{n+1}(t(x)) は n+1n+1 次で、最高次の係数は 2n(2/(b−a))n+12^n(2/(b-a))^{n+1}、零点はちょうど x0,…,xnx_0, \ldots, x_n です。ω−Tn+1(t(x))/(2n(2/(b−a))n+1)\omega - T_{n+1}(t(x))/\bigl(2^n(2/(b-a))^{n+1}\bigr) は Pn\mathcal{P}_n の元で n+1n+1 個の点 xjx_j で 00 なので、系 2.3 から ω(x)=2−n((b−a)/2)n+1Tn+1(t(x))\omega(x) = 2^{-n}((b-a)/2)^{n+1}T_{n+1}(t(x)) です。t(x)t(x) は [a,b][a,b] を [−1,1][-1,1] に写すので、補題 6.1 の 1・3 から max⁡[a,b]∣ω∣=2−n((b−a)/2)n+1=2((b−a)/4)n+1\max_{[a,b]}|\omega| = 2^{-n}((b-a)/2)^{n+1} = 2((b-a)/4)^{n+1} です。他の点の選び方でも ω\omega はモニックな n+1n+1 次式なので、定理 6.2 の 2 からこれ以上です。誤差の評価は、定理 4.1 の ξ\xi が [a,b][a,b] に入ることから従います。

補題 6.4(区間の外での Chebyshev 多項式の増大)✓

x≥1x \ge 1、w=x+x2−1w = x + \sqrt{x^2 - 1} とすると、k≥0k \ge 0 について

Tk(x)=wk+w−k2,wk2≤Tk(x)≤wkT_k(x) = \frac{w^k + w^{-k}}{2}, \qquad \frac{w^k}{2} \le T_k(x) \le w^k

である。また x≤−1x \le -1 では ∣Tk(x)∣=Tk(∣x∣)|T_k(x)| = T_k(|x|) である。

証明(補題 6.4)

w≥1w \ge 1 で w−1=x−x2−1w^{-1} = x - \sqrt{x^2 - 1} なので w+w−1=2xw + w^{-1} = 2x です。Sk=(wk+w−k)/2S_k = (w^k + w^{-k})/2 とおくと S0=1S_0 = 1、S1=xS_1 = x で、

2xSk−Sk−1=(w+w−1)(wk+w−k)−wk−1−w1−k2=wk+1+w−k−12=Sk+12xS_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}

なので、TkT_k と同じ漸化式と初期値を満たし、帰納法で Tk(x)=SkT_k(x) = S_k です。0<w−k≤1≤wk0 \lt w^{-k} \le 1 \le w^k から不等式が従います。最後の主張は 補題 6.1 の 2 の偶奇性です。

この補題は、区間 [−1,1][-1,1] の中では 11 以下に抑えられている TkT_k が、外に出ると wkw^k の速さで増えることを述べています。「区間の上では小さく、区間の外の 1 点では大きい」多項式が必要な場面(共役勾配法の収束の評価など)で、Chebyshev 多項式が最適な道具になる理由です(演習 7.4 を参照)。

定理 6.5(Chebyshev 点の Lebesgue 定数の上界)✓

n≥0n \ge 0 とし、[−1,1][-1,1] の Chebyshev 点 xj=cos⁡((2j+1)π/(2n+2))x_j = \cos\bigl((2j+1)\pi/(2n+2)\bigr)(j=0,…,nj = 0, \ldots, n)の Lebesgue 定数を Λn\Lambda_n とする。このとき

Λn≤2πlog⁡(2n+1)+1+4π\Lambda_n \le \frac{2}{\pi}\log(2n+1) + 1 + \frac{4}{\pi}

である。[a,b][a,b] の Chebyshev 点(系 6.3)でも同じ評価が成り立つ。

証明は長いので Appendix に置きます。最後の主張は、ℓj\ell_j が 1 次の変数変換 t(x)t(x) で [−1,1][-1,1] の場合の ℓj\ell_j に移ることから従います。等間隔点の指数的な増大と比べると、Chebyshev 点の Lebesgue 定数は n=1000n = 1000 でもこの評価から 7.27.2 以下です。

例 6.6(Lebesgue 定数を数値で求める)

格子(各小区間を 20 等分した点)の上で Lebesgue 関数の最大を求め、等間隔点では 命題 5.4 の下界と、Chebyshev 点では 定理 6.5 の上界と比べます。等間隔点の値は ℓj\ell_j の積に大きな数と小さな数が混ざるので、丸めの影響でないことを mpmath の 40 桁の計算で確かめます。格子上の最大は Λn\Lambda_n を下から近似する値です。

# 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}")
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=2n = 2 の場合は手で計算でき(演習 7.3)、等間隔点で 5/45/4、Chebyshev 点で 5/35/3 です。(b) 等間隔点の値は倍精度と 40 桁で表示の 5 桁まで一致し、丸め誤差の影響ではありません。(c) 等間隔点の値は下界 2n/(4n2)2^n/(4n^2) の上にあり、nn を 10 増やすごとに数百倍(20→3020 \to 30 で約 600600 倍、1 増やすごとにおよそ 1.91.9 倍)になって、下界の 2n2^n に近い速さで増えています。Chebyshev 点の値は上界の下にあり、nn を 10 増やしても 0.250.25〜0.410.41 ずつしか増えません。参考までに、よく引用される評価 (2/π)log⁡(n+1)+1(2/\pi)\log(n+1) + 1 は n=30n = 30 で 3.1863.186 で、表の 3.1493.149 はその下にもあります。

Lebesgue 定数と最良近似を組み合わせると、Chebyshev 点での補間が解析関数に対して幾何級数的に収束することが、導関数を使わずに示せます。

次の系の「Bernstein 楕円 EρE_\rho」は、複素数の平面で [−1,1][-1,1] を囲む、焦点が ±1\pm1 の楕円で、ρ\rho が大きいほど大きな楕円です。「正則」は複素数の変数で微分できることで、ff が EρE_\rho の内側でこの意味でなめらかなら、ρ\rho が大きいほど補間の誤差が速く減る、というのが系の内容です。

系 6.7(解析関数の Chebyshev 補間)✓

ρ>1\rho \gt 1 とし、ff は Bernstein 楕円 EρE_\rho で正則、[−1,1][-1,1] 上で実数値とする。pnp_n を [−1,1][-1,1] の Chebyshev 点での ff の補間多項式とすると、1<r<ρ1 \lt r \lt \rho と Mr=max⁡∂Er∣f∣M_r = \max_{\partial E_r}|f| について

max⁡[−1,1]∣f−pn∣≤(2+4π+2πlog⁡(2n+1))2Mr r−nr−1\max_{[-1,1]}|f - p_n| \le \Bigl(2 + \frac4\pi + \frac2\pi\log(2n+1)\Bigr)\frac{2M_r\,r^{-n}}{r - 1}

が成り立つ。

証明(系 6.7)

Chebyshev 展開の打ち切り(系 3.2[誤差保証つきの求積と関数近似]。係数の減衰は 定理 3.1[誤差保証つきの求積と関数近似])∑k=0nakTk\sum_{k=0}^{n}a_kT_k は Pn\mathcal{P}_n の元で、ff との差は [−1,1][-1,1] 上で 2Mrr−n/(r−1)2M_rr^{-n}/(r-1) 以下です。よって En(f)≤2Mrr−n/(r−1)E_n(f) \le 2M_rr^{-n}/(r-1) です。定理 5.2 の 2 と 定理 6.5 から主張を得ます。

対数の因子は幾何級数の減衰に比べて無視できるので、誤差はおよそ r−nr^{-n} の速さで減ります。rr は ρ\rho にいくらでも近く取れるので、収束の速さを決めるのは ff が正則でいられる最大の Bernstein 楕円の大きさ ρ\rho です。Runge の関数は ±i/5\pm i/5 に極をもち、ρ=1/5+1+1/25=(1+26)/5≈1.2198\rho = 1/5 + \sqrt{1 + 1/25} = (1 + \sqrt{26})/5 \approx 1.2198 です。

例 6.8(Chebyshev 補間の誤差の減り方)

exe^x には 系 6.3 の評価 e/(2n(n+1)!)e/(2^n(n+1)!) を、Runge の関数には 系 6.7 の予言する速さ ρ−n\rho^{-n} を当てはめます。

# 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}")
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) exe^x の誤差はすべて評価 e/(2n(n+1)!)e/(2^n(n+1)!) 以下で(比が 1 未満)、比がほぼ一定なので、誤差は評価式と同じ速さで減っています。比が 11 より小さいのは、評価が ∣f(n+1)∣|f^{(n+1)}| の最大値 ee を使っているためです。(b) Runge の関数の n=10,20,40n = 10, 20, 40 の誤差は、別の点の配置と並べた先の例(例 4.3)の Chebyshev 点の列と一致します。そこでは別の算法(chebinterpolate)とも一致していました。(c) 誤差に ρn\rho^n を掛けた値は 0.820.82 付近で一定になり、誤差がちょうど ρ−n\rho^{-n} の速さで減ること、すなわち nn を 10 増やすごとに ρ−10≈0.137\rho^{-10} \approx 0.137 倍になることを示しています。これは 系 6.7 の予言する速さ(任意の r<ρr \lt \rho について r−nr^{-n})と合っています。

注意 6.9(重心形の補間公式)

Lagrange 形を wj=1/∏k≠j(xj−xk)w_j = 1/\prod_{k \ne j}(x_j - x_k) を使って変形すると、pn(x)=∑jwjf(xj)x−xj/∑jwjx−xjp_n(x) = \sum_j \frac{w_j f(x_j)}{x - x_j} \Big/ \sum_j \frac{w_j}{x - x_j} と書けます(分母は定数 11 の補間から出ます)。これを重心形といい、重み wjw_j を一度計算すれば 1 点あたり O(n)O(n) の演算で値が求まり、点を足すときの更新も簡単です。Chebyshev 点では wjw_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)です。SciPy の BarycentricInterpolator はこの公式を使っています。

注意 6.10(最良近似と交代定理)

En(f)E_n(f) の下限は達成され(Pn\mathcal{P}_n は有限次元なので、最良近似多項式 q∗q^* が存在する)、q∗q^* はただ 1 つで、f−q∗f - q^* が [a,b][a,b] の n+2n+2 個の点で符号を交互に変えながら絶対値 En(f)E_n(f) をとることで特徴づけられます(Chebyshev の交代定理)。定理 6.2 はこの定理の f=xn+1f = x^{n+1} の場合にあたります。q∗q^* を数値で求める算法が Remez 法です。証明は T. J. Rivlin『An Introduction to the Approximation of Functions』(Dover, 1981)第 1 章、Trefethen『Approximation Theory and Approximation Practice』第 10 章にあります。本章では 定理 5.2 の 2 のように下限 En(f)E_n(f) を直接使うので、存在と一意性は必要ありません。

注意 6.11(スプライン)

次数を上げずに精度を上げるもう 1 つの方法は、区間を小区間に分けて各小区間で低次の多項式を使い、つなぎ目で滑らかさを課すことです。3 次スプライン(つなぎ目で C2C^2 級)は、刻み hh について誤差 O(h4)O(h^4) で、点を増やしても Runge 現象のような不安定さが起こりません。理論と算法は C. de Boor『A Practical Guide to Splines』(Springer, 1978、改訂版 2001)を参照してください。

演習 7.1易

f(x)=x3f(x) = x^3 の、点 0,1,2,30, 1, 2, 3 での差分商の表を定理 3.3の漸化式で作り、Newton 形が x3x^3 に一致すること、f[0,1,2,3]f[0,1,2,3] が定理 3.5の f′′′(ξ)/3!f'''(\xi)/3! に等しいことを確かめてください。

解答

値は 0,1,8,270, 1, 8, 27、1 階の差分商は 1,7,191, 7, 19、2 階は (7−1)/2=3(7-1)/2 = 3、(19−7)/2=6(19-7)/2 = 6、3 階は (6−3)/3=1(6-3)/3 = 1 です。Newton 形は

0+1⋅x+3x(x−1)+1⋅x(x−1)(x−2)=x+3x2−3x+x3−3x2+2x=x30 + 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′′′=6f''' = 6 なので f′′′(ξ)/3!=1=f[0,1,2,3]f'''(\xi)/3! = 1 = f[0,1,2,3] で、ξ\xi はどこでもかまいません。

演習 7.2標準

[0,1][0,1] の Chebyshev 点で exe^x を補間します。系 6.3 の評価で誤差 10−810^{-8} 以下が保証される最小の nn を求めてください。

解答

[0,1][0,1] 上 ∣f(n+1)∣≤e|f^{(n+1)}| \le e なので、評価は 2e(1/4)n+1/(n+1)!2e(1/4)^{n+1}/(n+1)! です。n=6n = 6 では 2e/(47⋅7!)=2e/82,575,360≈6.6×10−82e/(4^7 \cdot 7!) = 2e/82{,}575{,}360 \approx 6.6 \times 10^{-8} で足りず、n=7n = 7 では 2e/(48⋅8!)=2e/2,642,411,520≈2.1×10−92e/(4^8\cdot 8!) = 2e/2{,}642{,}411{,}520 \approx 2.1 \times 10^{-9} です。よって n=7n = 7(8 点)です。

演習 7.3標準

[−1,1][-1,1] の 3 点について Lebesgue 定数 Λ2\Lambda_2 を求めてください。(i) 等間隔点 −1,0,1-1, 0, 1。(ii) Chebyshev 点 −3/2,0,3/2-\sqrt3/2, 0, \sqrt3/2。どちらが小さいか、その理由も考えてください。

解答

(i) ℓ0=x(x−1)/2\ell_0 = x(x-1)/2、ℓ1=1−x2\ell_1 = 1 - x^2、ℓ2=x(x+1)/2\ell_2 = x(x+1)/2 です。対称性から [0,1][0,1] で考えると λ2=x(1−x)/2+(1−x2)+x(x+1)/2=1+x−x2\lambda_2 = x(1-x)/2 + (1 - x^2) + x(x+1)/2 = 1 + x - x^2 で、最大は x=1/2x = 1/2 での 5/45/4 です。

(ii) s=3/2s = \sqrt3/2 とすると ℓ0=x(x−s)/(2s2)\ell_0 = x(x-s)/(2s^2)、ℓ1=(s2−x2)/s2\ell_1 = (s^2 - x^2)/s^2、ℓ2=x(x+s)/(2s2)\ell_2 = x(x+s)/(2s^2) です。0≤x≤s0 \le x \le s では λ2=1+(xs−x2)/s2\lambda_2 = 1 + (xs - x^2)/s^2 で最大は x=s/2x = s/2 での 5/45/4、s≤x≤1s \le x \le 1 では ℓ1≤0\ell_1 \le 0 なので λ2=2x2/s2−1\lambda_2 = 2x^2/s^2 - 1 で、最大は x=1x = 1 での 8/3−1=5/38/3 - 1 = 5/3 です。よって Λ2=5/3\Lambda_2 = 5/3 です。

n=2n = 2 では等間隔点のほうが小さく、Chebyshev 点(第 1 種)は端点 ±1\pm1 を含まないので、端の近くで補間が外挿になって大きくなります。Chebyshev 点の利点は nn が大きいときの対数的な増え方にあり、数値の例(例 6.6)では n=5n = 5 ですでに逆転しています。

演習 7.4難

0<a<b0 \lt a \lt b、κ=b/a\kappa = b/a とします。p(0)=1p(0) = 1 を満たす kk 次以下の多項式 pp で

max⁡x∈[a,b]∣p(x)∣≤2(κ−1κ+1)k\max_{x \in [a,b]}|p(x)| \le 2\Bigl(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\Bigr)^k

となるものがあることを、補題 6.4 を使って示してください(この評価は共役勾配法の収束の評価の形です)。

解答

p(x)=Tk((b+a−2x)/(b−a))/Tk(x0)p(x) = T_k\bigl((b + a - 2x)/(b - a)\bigr)/T_k(x_0)、x0=(b+a)/(b−a)x_0 = (b+a)/(b-a) とおきます。x0>1x_0 \gt 1 なので 補題 6.4 から Tk(x0)>0T_k(x_0) \gt 0 で、p(0)=1p(0) = 1 です。x∈[a,b]x \in [a,b] なら (b+a−2x)/(b−a)∈[−1,1](b+a-2x)/(b-a) \in [-1,1] なので、補題 6.1 の 1 から分子の絶対値は 11 以下です。分母について、x0=(κ+1)/(κ−1)x_0 = (\kappa+1)/(\kappa-1)、x02−1=4κ/(κ−1)2x_0^2 - 1 = 4\kappa/(\kappa-1)^2 なので

w=x0+x02−1=κ+1+2κκ−1=(κ+1)2(κ−1)(κ+1)=κ+1κ−1w = 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}

で、補題 6.4 から Tk(x0)≥wk/2T_k(x_0) \ge w^k/2 です。よって max⁡[a,b]∣p∣≤2w−k\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

Appendix: Chebyshev 点の Lebesgue 定数の上界の証明

Section titled “Appendix: Chebyshev 点の Lebesgue 定数の上界の証明”

証明の流れは次のとおりです。まず ℓj\ell_j を θ\theta で表し、次に λn\lambda_n を余接(cot⁡\cot)の和で上から押さえます。最後に、xx に最も近い節点の項だけを別に扱い、残りの項を調和級数に似た和 ∑1/(2i−1)\sum 1/(2i-1) で評価して、それが対数で抑えられることを使います。

記号. N=n+1N = n + 1、θj=(2j+1)π/(2N)\theta_j = (2j+1)\pi/(2N)(j=0,…,nj = 0, \ldots, n)とし、x∈[−1,1]x \in [-1,1] を x=cos⁡θx = \cos\theta(θ∈[0,π]\theta \in [0,\pi])と書きます。証明に使う初等的な不等式を先に確かめます。平均値の定理(定理 3.3[平均値の定理とテイラーの定理])から、0<u<π/20 \lt u \lt \pi/2 で tan⁡u=u/cos⁡2c≥u\tan u = u/\cos^2 c \ge u(0<c<u0 \lt c \lt u)、したがって cot⁡u≤1/u\cot u \le 1/u です。cot⁡\cot は (0,π/2](0, \pi/2] で正、導関数 −1/sin⁡2u<0-1/\sin^2 u \lt 0 なので狭義減少です(系 3.4[平均値の定理とテイラーの定理])。y>0y \gt 0 で log⁡(1+y)=y/(1+c)≥y/(1+y)\log(1+y) = y/(1+c) \ge y/(1+y)(0<c<y0 \lt c \lt y)です。同じく平均値の定理から、実数 a,ba, b について ∣cos⁡Na−cos⁡Nb∣=N∣sin⁡Nc∣ ∣a−b∣≤N∣a−b∣|\cos Na - \cos Nb| = N|\sin Nc|\,|a - b| \le N|a - b|(cc は aa と bb の間)です。

証明(定理 6.5)

ℓj\ell_j の表示. 系 6.3 の証明のとおり ω=2−nTN\omega = 2^{-n}T_N で、積の微分から ω′(xj)=∏k≠j(xj−xk)\omega'(x_j) = \prod_{k \ne j}(x_j - x_k) なので、x≠xjx \ne x_j で ℓj(x)=ω(x)/((x−xj)ω′(xj))=TN(x)/((x−xj)TN′(xj))\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) です。補題 6.1 の 1 の両辺を θ\theta で微分すると TN′(cos⁡θ)sin⁡θ=Nsin⁡NθT_N'(\cos\theta)\sin\theta = N\sin N\theta で、sin⁡Nθj=sin⁡((2j+1)π/2)=(−1)j\sin N\theta_j = \sin\bigl((2j+1)\pi/2\bigr) = (-1)^j から ∣TN′(xj)∣=N/sin⁡θj|T_N'(x_j)| = N/\sin\theta_j です。よって x=cos⁡θ≠xjx = \cos\theta \ne x_j で

∣ℓj(x)∣=∣cos⁡Nθ∣N⋅sin⁡θj∣cos⁡θ−cos⁡θj∣|\ell_j(x)| = \frac{|\cos N\theta|}{N}\cdot\frac{\sin\theta_j}{|\cos\theta - \cos\theta_j|}

です。xx がどれかの節点なら λn(x)=1\lambda_n(x) = 1 なので、以下 xx は節点でないとします。

余接への分解. α=(θj+θ)/2\alpha = (\theta_j + \theta)/2、β=(θj−θ)/2\beta = (\theta_j - \theta)/2 とおくと、cos⁡θ−cos⁡θj=2sin⁡αsin⁡β\cos\theta - \cos\theta_j = 2\sin\alpha\sin\beta、sin⁡θj=sin⁡(α+β)=sin⁡αcos⁡β+cos⁡αsin⁡β\sin\theta_j = \sin(\alpha + \beta) = \sin\alpha\cos\beta + \cos\alpha\sin\beta なので

sin⁡θjcos⁡θ−cos⁡θj=12(cot⁡β+cot⁡α)\frac{\sin\theta_j}{\cos\theta - \cos\theta_j} = \frac12\bigl(\cot\beta + \cot\alpha\bigr)

です(0<α<π0 \lt \alpha \lt \pi、β∉πZ\beta \notin \pi\mathbb{Z})。∣cot⁡∣|\cot| は偶関数で周期 π\pi なので ∣cot⁡α∣=∣cot⁡((θ−(2π−θj))/2)∣|\cot\alpha| = |\cot\bigl((\theta - (2\pi - \theta_j))/2\bigr)| です。{θj}∪{2π−θj}\{\theta_j\} \cup \{2\pi - \theta_j\} はちょうど sm=(2m+1)π/(2N)s_m = (2m+1)\pi/(2N)(m=0,…,2N−1m = 0, \ldots, 2N-1)の全体なので、vm=(θ−sm)/2v_m = (\theta - s_m)/2 とおくと

λn(x)≤∣cos⁡Nθ∣2N∑m=02N−1∣cot⁡vm∣\lambda_n(x) \le \frac{|\cos N\theta|}{2N}\sum_{m=0}^{2N-1}|\cot v_m|

です。

点の並び. 実数 vv と πZ\pi\mathbb{Z} の距離を δ(v)∈[0,π/2]\delta(v) \in [0, \pi/2] と書くと、∣cot⁡v∣=cot⁡δ(v)|\cot v| = \cot\delta(v) です。vmv_m は公差 −h-h(h=π/(2N)h = \pi/(2N)、2Nh=π2Nh = \pi)の 2N2N 項の等差数列なので、π\pi を法として周長 π\pi の円周上に間隔 hh で並ぶ 2N2N 点です。δ(vm0)=d\delta(v_{m_0}) = d が最小となる m0m_0 を取ると d≤h/2d \le h/2 です。xx は節点でないので d>0d \gt 0 です。m≠m0m \ne m_0 では vm−vm0=khv_m - v_{m_0} = kh で、kk を 2N2N で割った余り k′k' は 1,…,2N−11, \ldots, 2N-1 を 1 回ずつとり、δ(kh)=min⁡(k′,2N−k′) h\delta(kh) = \min(k', 2N - k')\,h です。vm0=pπ+ev_{m_0} = p\pi + e(p∈Zp \in \mathbb{Z}、∣e∣=d|e| = d)と書くと δ(vm)=δ(e+kh)≥δ(kh)−d\delta(v_m) = \delta(e + kh) \ge \delta(kh) - d なので、i=min⁡(k′,2N−k′)i = \min(k', 2N - k') として δ(vm)≥(i−1/2)h\delta(v_m) \ge (i - 1/2)h です。ii は 1,…,N−11, \ldots, N-1 を 2 回ずつ、NN を 1 回とります。cot⁡\cot の減少性から

∑m≠m0∣cot⁡vm∣≤2∑i=1Ncot⁡(2i−1)π4N≤2∑i=1N4N(2i−1)π\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⁡(Nsm0)=cos⁡((2m0+1)π/2)=0\cos(Ns_{m_0}) = \cos\bigl((2m_0+1)\pi/2\bigr) = 0 で、θ−(sm0+2pπ)=2e\theta - (s_{m_0} + 2p\pi) = 2e なので、∣cos⁡Nθ∣=∣cos⁡Nθ−cos⁡N(sm0+2pπ)∣≤2Nd|\cos N\theta| = |\cos N\theta - \cos N(s_{m_0} + 2p\pi)| \le 2Nd です。よって

∣cos⁡Nθ∣2N∣cot⁡vm0∣≤2Nd2Ncot⁡d≤d⋅1d=1\frac{|\cos N\theta|}{2N}|\cot v_{m_0}| \le \frac{2Nd}{2N}\cot d \le d \cdot \frac1d = 1

です。

まとめ. 残りの項は ∣cos⁡Nθ∣≤1|\cos N\theta| \le 1 と上の評価から 12N⋅8Nπ∑i=1N12i−1=4π∑i=1N12i−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≥2i \ge 2 では y=2/(2i−3)y = 2/(2i-3) として log⁡2i−12i−3=log⁡(1+y)≥y1+y=22i−1\log\frac{2i-1}{2i-3} = \log(1+y) \ge \frac{y}{1+y} = \frac{2}{2i-1} なので、

∑i=1N12i−1≤1+12∑i=2Nlog⁡2i−12i−3=1+12log⁡(2N−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+12N - 1 = 2n+1 なので λn(x)≤1+4π+2πlog⁡(2n+1)\lambda_n(x) \le 1 + \frac4\pi + \frac2\pi\log(2n+1) で、xx について最大をとれば主張を得ます。

[a,b][a,b] の Chebyshev 点 xj=a+b2+b−a2x^jx_j = \frac{a+b}{2} + \frac{b-a}{2}\hat x_j では、t(x)=(2x−a−b)/(b−a)t(x) = (2x - a - b)/(b - a) について ℓj(x)=ℓ^j(t(x))\ell_j(x) = \hat\ell_j(t(x))(ℓ^j\hat\ell_j は [−1,1][-1,1] の点 x^j\hat x_j の Lagrange 基底)です。x−xk=b−a2(t(x)−x^k)x - x_k = \frac{b-a}{2}(t(x) - \hat x_k) の比をとれば因子 b−a2\frac{b-a}{2} が約分されるからです。tt は [a,b][a,b] を [−1,1][-1,1] の上へ写すので、Lebesgue 定数は等しくなります。

このテーマの 100 本ノック(50 問)

執筆: Claude Code(Opus 5.5)、Claude Code(Sonnet 5.5)・監査 2 回(Claude Code)・最終 2026-10-12

この記事の誤りを報告する ・運営: 夢現技研合同会社 ・料金 ・利用条件 ・特定商取引法に基づく表記

© 2026 夢現技研合同会社 ・コード例は MIT ライセンスです。