誤差保証つきの求積と関数近似：Chebyshev 展開・Bernstein 楕円・区間 Taylor 法

0. この記事の要点

$e^{-x^2}$ の積分には初等関数で書ける原始関数がなく、値は数値計算で求めるしかありません。ふつうの数値計算は「だいたいこの値」を返すだけですが、定理の証明に使うには「真の値は必ずこの区間の中にある」と言い切れる必要があります。この記事は、積分の値を、誤差まで込みで確実に閉じ込める 2 つの方法を作ります。読み終えると、必要な点数の見積もりと、積分を含む区間の計算ができるようになります。

• $[-1,1]$ の近くで正則な関数は Chebyshev 多項式で展開でき、係数は $|a_k|\le2M_rr^{-k}$ と幾何級数的に減衰します。$r$ は関数が正則でいられる Bernstein 楕円の大きさ、$M_r$ はその楕円の上での $|f|$ の最大値です。
• 証明の鍵は Joukowski 変換 $x=\tfrac12(w+w^{-1})$ です。これで Chebyshev 展開は円環上の Laurent 展開になり、係数の評価は Cauchy の評価そのものになります。
• $n$ 点 Gauss–Legendre 求積は $2n-1$ 次まで厳密なので、誤差は $2n$ 次以上の Chebyshev 係数だけで決まり、$r^{-2n}$ の速さで小さくなります。上界は $M_r$ を区間演算で評価すれば計算できます。
• 区間 Taylor 法は、Taylor 係数を漸化式で計算し、剰余項を区間全体で評価した係数で包みます。四則と単調性だけで $\int_0^1e^{-x^2}dx$ を幅 $10^{-14}$ 程度の区間に閉じ込められます。

この記事で使う事前知識

• 区間演算: 実数の代わりに $[\text{下端},\text{上端}]$ の形の区間で計算し、結果が真の値を必ず含むようにする方法。→ 浮動小数と区間演算
• 丸め: 実数を計算機で表せる数（浮動小数）に置き換える操作。区間演算では端点を外向きに丸めます。→ 丸め
• Laurent 展開: 円環上で正則な関数を、負べきも含む $\sum c_kw^k$ の形に展開すること。→ Laurent の定理
• Cauchy の評価: 展開の係数の大きさを、円周上の $|f|$ の最大値と半径で抑える不等式。→ コーシーの評価式
• Gauss 求積: $n$ 点で $2n-1$ 次以下の多項式を厳密に積分する公式。→ Gauss 求積
• 整関数: 複素平面全体で正則な関数。→ 整関数

1. 動機

積分の値が定理の一部になる場面は多くあります。固有値の評価で行列成分を積分で定めるとき、特殊関数の値で不等式を確かめるときなどです。数値積分の公式には誤差項がありますが、教科書にある $f^{(2n)}(\xi)$ の形の誤差項は、高階導関数の上界がわからないと使えません。

本章では、計算機で実際に評価できる上界を 2 通り作ります。1 つは関数を複素平面に延ばし、楕円の上での $|f|$ の最大値だけで誤差を抑える方法です（第 3・4 節）。高階導関数は要りません。もう 1 つは Taylor 係数を漸化式で計算し、その剰余を前章の区間演算で包む方法です（第 5 節）。前者はどれだけの点数で足りるかの設計に、後者は値を含む区間そのものの計算に向いています。

2. Chebyshev 多項式と Joukowski 変換

定義 2.1（Chebyshev 多項式と Bernstein 楕円）↗

$T_0(x)=1$, $T_1(x)=x$, $T_{k+1}(x)=2xT_k(x)-T_{k-1}(x)$ で定まる多項式 $T_k$ を Chebyshev 多項式という。$J(w)=\tfrac12(w+w^{-1})$（$w\ne0$）を Joukowski 変換という。$r\gt1$ に対し、焦点 $\pm1$、長半径 $\tfrac12(r+r^{-1})$、短半径 $\tfrac12(r-r^{-1})$ の楕円 $\partial E_r$ と、その内部の開領域 $E_r$ を Bernstein 楕円という。

次の補題のねらいは、$w$ 平面の円を $J$ で写すと、$x$ 平面では $[-1,1]$ を焦点とする楕円や線分になる、と見ることです。半径 $r$ の円と半径 $r^{-1}$ の円は $J(w)=J(w^{-1})$ のため同じ楕円に重なり、単位円は $[-1,1]$ を往復します。この見方で Chebyshev 展開を円の上の展開として扱えるようになります。

補題 2.2（Joukowski 変換の性質）✓

$\rho\gt1$ とする。

1. すべての $k\ge0$ と $w\ne0$ について $T_k(J(w))=\tfrac12(w^k+w^{-k})$。とくに $T_k(\cos\theta)=\cos k\theta$ で、$[-1,1]$ 上 $|T_k|\le1$。
2. $1\lt r\lt\rho$ のとき、$J$ は円 $|w|=r$ と円 $|w|=r^{-1}$ をともに楕円 $\partial E_r$ の上へ写し、単位円を $[-1,1]$ の上へ写す。
3. $J$ は円環 $A_\rho=\{\rho^{-1}\lt|w|\lt\rho\}$ を $E_\rho$ の中へ写す。

証明（補題 2.2）

1. $k=0,1$ では定義どおりです。$k$, $k-1$ で成り立てば

$$2J(w)\,T_k(J(w))-T_{k-1}(J(w))=\tfrac12(w+w^{-1})(w^k+w^{-k})-\tfrac12(w^{k-1}+w^{1-k})=\tfrac12(w^{k+1}+w^{-k-1})$$

です。$w=e^{i\theta}$ とおくと $J(w)=\cos\theta$ で、右辺は $\cos k\theta$ になります。

2. $w=se^{i\theta}$ とおくと $J(w)=\tfrac12(s+s^{-1})\cos\theta+\tfrac i2(s-s^{-1})\sin\theta$ です。$s=r$ なら長半径 $a_r=\tfrac12(r+r^{-1})$、短半径 $b_r=\tfrac12(r-r^{-1})$ の楕円を 1 周し、$a_r^2-b_r^2=1$ なので焦点は $\pm1$ です。$s=r^{-1}$ では $J(w)=J(w^{-1})$ から同じ楕円、$s=1$ では $[-1,1]$ です。

3. $a_r$, $b_r$ は $r\gt1$ で狭義増加です。$1\lt s\lt\rho$ のとき $J(w)=x+iy$ は $x^2/a_s^2+y^2/b_s^2=1$ を満たすので、$x^2/a_\rho^2+y^2/b_\rho^2\lt1$、すなわち $J(w)\in E_\rho$ です。$s=1$ なら $J(w)\in[-1,1]\subset E_\rho$、$s\lt1$ なら $J(w)=J(w^{-1})$ に帰着します。

∎

例 6.4［直交多項式］ にあるとおり、$T_k$ は重み $1/\sqrt{1-x^2}$ に関する直交多項式でもあります。本章では直交性は使わず、補題 2.2 の 1 だけを使います。

3. Chebyshev 展開の係数の減衰

定理 3.1（Bernstein 楕円と Chebyshev 係数）✓

$\rho\gt1$ とし、$f$ は $E_\rho$ で正則で $[-1,1]$ 上で実数値とする。このとき実数 $a_0,a_1,\ldots$ があって

$$f(x)=\sum_{k=0}^\infty a_kT_k(x)\qquad(x\in[-1,1])$$

が絶対かつ一様に収束し、任意の $1\lt r\lt\rho$ について、$M_r=\max_{z\in\partial E_r}|f(z)|$ とおくと

$$|a_0|\le M_r,\qquad|a_k|\le2M_rr^{-k}\quad(k\ge1)$$

が成り立つ。

証明の方針は次のとおりです。$f(J(w))$ を $w$ の Laurent 展開にし、係数を円周上の積分で書きます。半径 $r$ の円では $k\ge0$ の係数が $r^{-k}$ で、半径 $r^{-1}$ の円では $k\lt0$ の係数が同じ形で抑えられます。2 つの円は同じ楕円 $\partial E_r$ に写るので、どちらも同じ $M_r$ が使えます。

証明（定理 3.1）

$F(w)=f(J(w))$ とおくと、補題 2.2 の 3 から $F$ は円環 $A_\rho$ で正則です。Laurent 展開（定理 6.1［正則関数の強力な性質］）により、$A_\rho$ 上で $F(w)=\sum_{k\in\mathbb{Z}}c_kw^k$ と書け、係数は $A_\rho$ 内の任意の円 $|w|=s$ 上の積分

$$c_k=\frac1{2\pi i}\int_{|w|=s}\frac{F(w)}{w^{k+1}}\,dw$$

で与えられます。

評価。 $s=r$ とすると、補題 2.2 の 2 から円 $|w|=r$ 上で $|F(w)|\le M_r$ なので、$|c_k|\le\frac1{2\pi}\cdot2\pi r\cdot M_rr^{-k-1}=M_rr^{-k}$ です。$s=r^{-1}$ とすると同様に $|c_k|\le M_rr^{k}$ です。よって $|c_k|\le M_rr^{-|k|}$ がすべての $k\in\mathbb{Z}$ で成り立ちます。

対称性。 $s=1$ とし $w=e^{i\theta}$ と置換すると $c_k=\frac1{2\pi}\int_{-\pi}^{\pi}f(\cos\theta)e^{-ik\theta}\,d\theta$ です。$f(\cos\theta)$ は $\theta$ の偶関数で実数値なので、$c_k=\frac1{2\pi}\int_{-\pi}^{\pi}f(\cos\theta)\cos k\theta\,d\theta$ は実数で $c_{-k}=c_k$ です。

展開。 単位円上で $\sum_k|c_k|\le M_r+2M_r\sum_{k\ge1}r^{-k}\lt\infty$ なので、$F(e^{i\theta})=c_0+\sum_{k\ge1}c_k(e^{ik\theta}+e^{-ik\theta})$ は絶対一様収束します。$a_0=c_0$, $a_k=2c_k$（$k\ge1$）とおけば、補題 2.2 の 1 により $f(\cos\theta)=\sum_ka_kT_k(\cos\theta)$ で、$|a_k|=2|c_k|\le2M_rr^{-k}$ です。

∎

系 3.2（打ち切り誤差）✓

定理 3.1 の仮定のもとで、$x\in[-1,1]$ について

$$\Bigl|f(x)-\sum_{k=0}^na_kT_k(x)\Bigr|\le\frac{2M_rr^{-n}}{r-1}.$$

証明（系 3.2）

$|T_k|\le1$ から、左辺は $\sum_{k\gt n}|a_k|\le2M_r\sum_{k\gt n}r^{-k}=2M_rr^{-n}/(r-1)$ 以下です。

∎

この評価は、高階導関数をまったく使いません。必要なのは楕円 $\partial E_r$ の上での $|f|$ の上界で、これは楕円を小さな弧に分けて複素数の区間演算で評価すれば得られます。$r$ を大きくとると $r^{-k}$ は速く減りますが $M_r$ は増えるので、両者の釣り合いで最良の $r$ が決まります。

Bernstein 楕円 ∂E_r（r = 1.5, 2）と、Runge の関数 1/(1+25x²) の極 ±i/5 を通る楕円（破線、r ≈ 1.22）。f が正則でいられる楕円が大きいほど、Chebyshev 係数は速く減衰する。

楕円は $r$ が大きいほど横にも縦にも広がります。$f$ が正則でいられるのは極を含まない楕円の内側だけなので、虚軸上の極 $\pm i/5$ に楕円がぶつかる大きさが、使える $\rho$ の上限です。次の例では、その大きさを短半径から求めています。

例 3.3（Runge の関数と Gauss 関数）

$f(x)=1/(1+25x^2)$ は $\pm i/5$ に極をもちます。$E_\rho$ がこの極を含まない最大の $\rho$ は、短半径 $\tfrac12(\rho-\rho^{-1})=\tfrac15$ から $\rho=\tfrac15+\sqrt{1+\tfrac1{25}}\approx1.2198$ です。係数は $1.22^{-k}$ 程度の遅さでしか減らず、$10^{-10}$ の精度に $k\approx120$ 項が要ります。一方 $f(x)=e^{-x^2}$ は整関数で、$z=a_r\cos\theta+ib_r\sin\theta$ の上で $|e^{-z^2}|=e^{-\mathrm{Re}\,z^2}=e^{-a_r^2\cos^2\theta+b_r^2\sin^2\theta}\le e^{b_r^2}$ です。どの $r$ でも評価が使えるので、係数は幾何級数よりも速く減衰します。

4. Gauss–Legendre 求積の誤差

$[-1,1]$ 上の Lebesgue 測度に対する直交多項式（Legendre 多項式）の零点 $x_1\lt\cdots\lt x_n$ と重み $w_i$ による $G_n(f)=\sum_iw_if(x_i)$ を $n$ 点 Gauss–Legendre 求積といいます。定理 6.1［直交多項式］ により、$G_n$ は $2n-1$ 次以下の多項式に対して厳密で、重みは正です。$I(f)=\int_{-1}^1f(x)\,dx$ と書きます。

補題 4.1（節点と重みの対称性）✓

$x_{n+1-i}=-x_i$, $w_{n+1-i}=w_i$ であり、$\sum_iw_i=2$ である。したがって奇関数 $g$ について $G_n(g)=0=I(g)$。

証明（補題 4.1）

$P_n$ を $n$ 次のモニックな直交多項式とすると、$q(x)=(-1)^nP_n(-x)$ もモニックで、置換 $x\mapsto-x$ が測度 $dx$ を保つので、$n-1$ 次以下の多項式すべてと直交します。直交多項式の一意性（定理 2.2［直交多項式］）から $q=P_n$ で、零点の集合は $x\mapsto-x$ で保たれます。Lagrange 基底は $\ell_{n+1-i}(x)=\ell_i(-x)$ を満たすので、$w_{n+1-i}=\int\ell_i(-x)\,dx=w_i$ です。$\sum_iw_i=G_n(1)=I(1)=2$ は厳密性からです。最後の主張は、奇関数の値が対称な節点で打ち消し合うことからです。

∎

補題 4.2（Chebyshev 多項式の積分）✓

$k$ が奇数なら $I(T_k)=0$、偶数なら $I(T_k)=\dfrac{2}{1-k^2}$ である。

証明（補題 4.2）

漸化式から $T_k$ は $k$ と同じ偶奇の多項式なので、奇数では $0$ です。偶数 $k$ では $x=\cos\theta$ と置換して

$$I(T_k)=\int_0^\pi\cos k\theta\sin\theta\,d\theta=\frac12\int_0^\pi\bigl(\sin(k+1)\theta-\sin(k-1)\theta\bigr)d\theta=\frac12\Bigl(\frac2{k+1}-\frac2{k-1}\Bigr)=\frac2{1-k^2}$$

です（$k\pm1$ は奇数なので $\int_0^\pi\sin m\theta\,d\theta=2/m$。$k=0$ でも成り立ちます）。

∎

定理 4.3（Gauss–Legendre 求積の誤差の上界）✓

定理 3.1 の仮定のもとで、$1\lt r\lt\rho$ と $n\ge1$ について

$$|I(f)-G_n(f)|\le\frac{16n^2}{4n^2-1}\cdot\frac{M_r\,r^{2-2n}}{r^2-1}.$$

証明の方針は、$f$ を Chebyshev 展開して項ごとに誤差を見ることです。$2n-1$ 次以下の項は求積が厳密なので誤差が出ず、奇数次の項は対称性で消えます。残る偶数次で $2n$ 次以上の項だけを、係数の減衰で足し上げます。

証明（定理 4.3）

$E=I-G_n$ とおきます。$f=\sum_ka_kT_k$ は $[-1,1]$ で一様収束するので、項別積分（定理 5.2［関数列と一様収束］）により $I(f)=\sum_ka_kI(T_k)$ です。$G_n$ は有限個の点の値の和なので $G_n(f)=\sum_ka_kG_n(T_k)$ です。よって $E(f)=\sum_ka_kE(T_k)$ です。

$k\le2n-1$ なら厳密性から $E(T_k)=0$、$k$ が奇数なら 補題 4.1 から $E(T_k)=0$ です。残るのは偶数 $k\ge2n$ で、補題 4.2、$|T_k|\le1$、重みが正で和が $2$ であることから

$$|E(T_k)|\le|I(T_k)|+\sum_iw_i|T_k(x_i)|\le\frac2{k^2-1}+2\le\frac2{4n^2-1}+2=\frac{8n^2}{4n^2-1}$$

です。定理 3.1 の $|a_k|\le2M_rr^{-k}$ と合わせて

$$|E(f)|\le\frac{8n^2}{4n^2-1}\cdot2M_r\sum_{j\ge n}r^{-2j}=\frac{16n^2}{4n^2-1}\cdot\frac{M_rr^{-2n}}{1-r^{-2}}$$

で、これは主張の右辺に等しいです。

∎

誤差が $r^{-2n}$ で減ること、つまり点数を 1 増やすごとに $r^2$ 倍よくなることが Gauss 求積の強さです。同じ点数の台形則や Simpson 則は多項式の速さでしか減りません（周期関数を除く）。

例 4.4（誤差の上界と実際の誤差）

$f(x)=e^{-x^2}$ では 例 3.3 から $M_r\le e^{b_r^2}$ です。各 $n$ で $r$ を数値的に動かして上界を最小にし、NumPy の leggauss で計算した実際の誤差と比べると次のようになります（真値は $\sqrt\pi\,\mathrm{erf}(1)$）。

$n$	実際の誤差	定理 4.3 の上界
2	$6.1\times10^{-2}$	$3.5\times10^{-1}$
4	$3.1\times10^{-4}$	$2.2\times10^{-3}$
6	$6.5\times10^{-7}$	$5.4\times10^{-6}$
8	$7.3\times10^{-10}$	$6.9\times10^{-9}$
10	$5.0\times10^{-13}$	$5.3\times10^{-12}$

上界は実際の誤差の 10 倍程度で、減り方の速さはよく一致しています。

注意 4.5

上の表の「実際の誤差」は浮動小数で計算したもので、証明ではありません。求積の和そのものを区間で計算するには、節点と重みも区間で包む必要があります。節点は $P_n$ の零点なので第 1 章の区間 Newton 法で包め、重みは 定理 6.1［直交多項式］ の 2 の公式を区間で評価すれば包めます。定理 4.3 の上界も、$M_r$ を区間演算で上から評価し、残りを上向き丸めで計算して初めて厳密になります。

5. 区間 Taylor 法

積分を含む区間を直接計算する方法に移ります。関数 $f$ の点 $c$ での Taylor 係数を $(f)_k=f^{(k)}(c)/k!$ と書きます。

補題 5.1（Taylor 係数の漸化式）✓

$f$, $g$ が $c$ の近傍で $k$ 回微分可能なら $(fg)_k=\sum_{j=0}^k(f)_j(g)_{k-j}$ である。また $g=e^f$ なら $(g)_0=e^{(f)_0}$ で、$k\ge1$ について

$$(g)_k=\frac1k\sum_{j=1}^kj\,(f)_j\,(g)_{k-j}.$$

証明（補題 5.1）

Leibniz の公式 $(fg)^{(k)}=\sum_j\binom kjf^{(j)}g^{(k-j)}$（積の微分法と $\binom kj+\binom k{j-1}=\binom{k+1}j$ による $k$ の帰納法）を $k!$ で割れば第 1 式です。$g=e^f$ なら $g'=f'g$ です。$(g')_{k-1}=g^{(k)}(c)/(k-1)!=k(g)_k$、同様に $(f')_{i}=(i+1)(f)_{i+1}$ なので、第 1 式を $f'g$ の $k-1$ 次の係数に使うと

$$k(g)_k=\sum_{i=0}^{k-1}(i+1)(f)_{i+1}(g)_{k-1-i}=\sum_{j=1}^kj\,(f)_j(g)_{k-j}$$

です。

∎

この漸化式は四則と $\exp$ だけでできた式です。したがって 定理 4.4［浮動小数と区間演算］ により、$c$ の代わりに区間 $C$ を入れて区間で計算すれば、すべての $c\in C$ についての $(f)_k$ を含む区間が得られます。ただし $\exp$ には、計算機の exp ではなく、証明つきの包含が必要です（第 1 章の注意）。

補題 5.2（指数関数の包含）✓

$x\le0$ と $K\ge0$ について $\Bigl|e^x-\sum_{k=0}^K\dfrac{x^k}{k!}\Bigr|\le\dfrac{|x|^{K+1}}{(K+1)!}$ である。

証明（補題 5.2）

$x=0$ なら両辺 $0$ です。$x\lt0$ なら Taylor の定理（定理 5.3［平均値の定理とテイラーの定理］）により、剰余は $e^\xi x^{K+1}/(K+1)!$（$x\lt\xi\lt0$）で、$0\lt e^\xi\lt1$ です。

∎

区間 $[p,q]\subset(-\infty,0]$ に対しては、$\exp$ の単調性から $e^{[p,q]}\subset[\,e^p\text{ の包含の下端},\ e^q\text{ の包含の上端}\,]$ とします。

次の定理の考え方は、$f$ を $c$ のまわりの $n-1$ 次の Taylor 多項式と剰余に分けることです。多項式の部分は中心での係数 $G_k$ で、剰余は $[L,U]$ で、それぞれ積分を区間として包みます。$n$ を偶数にするのは、剰余に掛かる $(x-c)^n$ が負にならず、$L$ と $U$ で挟む向きが変わらないようにするためです。

定理 5.3（区間 Taylor 法による積分の包含）✓

$n\ge2$ を偶数、$X=[c-h,c+h]$（$h\gt0$）とし、$f$ は $X$ を含む開区間で $C^n$ 級とする。区間 $G_k$ が $(f)_k=f^{(k)}(c)/k!$ を含み（$k=0,1,\ldots,n-1$）、区間 $[L,U]$ が $\{f^{(n)}(\xi)/n!:\xi\in X\}$ を含むとする。このとき

$$\int_{c-h}^{c+h}f(x)\,dx\ \in\ \sum_{\substack{0\le k\lt n\\ k\text{ は偶数}}}G_k\cdot\frac{2h^{k+1}}{k+1}+[L,U]\cdot\frac{2h^{n+1}}{n+1}.$$

証明（定理 5.3）

$r(x)=f(x)-\sum_{k\lt n}(f)_k(x-c)^k$ とおきます。$r$ は連続です。$x\ne c$ なら Taylor の定理（定理 5.3［平均値の定理とテイラーの定理］、$n-1$ 次まで展開）により $r(x)=\frac{f^{(n)}(\xi_x)}{n!}(x-c)^n$ となる $\xi_x\in X$ があり、$n$ が偶数なので $(x-c)^n\ge0$ です。よって

$$L(x-c)^n\le r(x)\le U(x-c)^n\qquad(x\in X)$$

で（$x=c$ では 3 辺とも $0$）、積分の単調性（命題 4.1［積分の基本定理と定積分］ の (2)）から $\int_Xr\in[L,U]\cdot\frac{2h^{n+1}}{n+1}$ です。また $\int_X(x-c)^kdx=\frac{h^{k+1}-(-h)^{k+1}}{k+1}$ は $k$ が奇数なら $0$、偶数なら $\frac{2h^{k+1}}{k+1}$ です。$\int_Xf=\sum_{k\lt n}(f)_k\int_X(x-c)^kdx+\int_Xr$ の各項がそれぞれの区間に入るので、区間演算の定義から和も右辺に入ります。

∎

剰余の区間 $[L,U]$ は、補題 5.1 の漸化式を $C=X$ として区間で評価した $n$ 次の係数で取れます。区間 $X$ 全体での評価なので過大評価されますが、$h^{n+1}$ が掛かるので、小区間に分ければ影響は小さくなります。$[a,b]$ 全体の積分は、小区間ごとの包含の和です。

次の例のコードは、(1) nextafter（浮動小数で次に大きい、または小さい数を返す関数）で端点を 1 つずつ外へ広げる区間の型 I、(2) 証明つきの $e^x$ の包含 exp_point、(3) 漸化式による Taylor 係数 taylor_coeffs、(4) 小区間ごとに定理を使って足し合わせる integrate01、の 4 つからなります。細部が追えなくても、4 を読めば定理の式がそのまま実装されていることがわかります。

例 5.4（∫₀¹ exp(−x²) dx の包含）

$f=e^{g}$, $g(x)=-x^2$ とすると、$g$ の Taylor 係数は $(g)_0=-c^2$, $(g)_1=-2c$, $(g)_2=-1$、それ以外は $0$ です。$[0,1]$ を $m$ 等分（$m$ は 2 の冪にして、端点・中点・半幅を浮動小数で厳密に表す）し、各小区間に 定理 5.3 を使います。区間の型 I は第 1 章のもの（四則を nextafter で外へ広げる）に、符号反転と平方を足したものです。

import math

INF = math.inf
def down(x): return math.nextafter(x, -INF)
def up(x):   return math.nextafter(x, INF)

class I:
    """第 1 章の区間（最近点丸めの結果を nextafter で 1 つずつ外へ広げる）"""
    def __init__(self, lo, hi=None):
        self.lo, self.hi = lo, (lo if hi is None else hi)
        assert self.lo <= self.hi
    def __add__(self, o): return I(down(self.lo + o.lo), up(self.hi + o.hi))
    def __sub__(self, o): return I(down(self.lo - o.hi), up(self.hi - o.lo))
    def __neg__(self):    return I(-self.hi, -self.lo)
    def __mul__(self, o):
        p = [self.lo * o.lo, self.lo * o.hi, self.hi * o.lo, self.hi * o.hi]
        return I(down(min(p)), up(max(p)))
    def __truediv__(self, o):
        assert not (o.lo <= 0 <= o.hi)
        q = [self.lo / o.lo, self.lo / o.hi, self.hi / o.lo, self.hi / o.hi]
        return I(down(min(q)), up(max(q)))
    def sqr(self):                                   # {x^2 : x in X}
        a, b = sorted([abs(self.lo), abs(self.hi)])
        lo = 0.0 if self.lo <= 0 <= self.hi else max(0.0, down(a * a))
        return I(lo, up(b * b))

def exp_point(x, K=25):
    """浮動小数 x <= 0 について e^x を含む区間。math.exp は使わない。"""
    X, term, s = I(x), I(1.0), I(1.0)
    for k in range(1, K + 1):
        term = term * X / I(float(k))                # x^k / k!
        s = s + term
    rem = term * X / I(float(K + 1))                 # x^{K+1} / (K+1)!
    r = max(abs(rem.lo), abs(rem.hi))
    return s + I(-r, r)                              # 剰余の上界（e^ξ <= 1）

def exp_interval(X):                                 # X ⊂ (-inf, 0]、exp は単調増加
    return I(exp_point(X.lo).lo, exp_point(X.hi).hi)

def taylor_coeffs(C, n):
    """f(x) = exp(-x^2) の Taylor 係数 f^(k)(c)/k! (c ∈ C) を含む区間、k = 0..n。"""
    g = [-C.sqr(), -(I(2.0) * C), I(-1.0)]           # g = -x^2 の係数（3 次以上は 0）
    f = [exp_interval(g[0])]
    for k in range(1, n + 1):                        # k f_k = Σ_j j g_j f_{k-j}
        s = I(0.0)
        for j in range(1, min(k, 2) + 1):
            s = s + I(float(j)) * g[j] * f[k - j]
        f.append(s / I(float(k)))
    return f

def integrate01(m, n):
    """∫_0^1 exp(-x^2) dx を含む区間。m は 2 の冪、n は偶数。"""
    total = I(0.0)
    h = 1.0 / (2 * m)                                # 小区間の半幅（厳密）
    for i in range(m):
        c = (2 * i + 1) * h                          # 中点（厳密）
        fc = taylor_coeffs(I(c), n)                  # 中点での係数
        fX = taylor_coeffs(I(c - h, c + h), n)       # 小区間全体での係数（剰余用）
        w = 2.0 * h                                  # 2h^{k+1}（2 の冪なので積は厳密）
        for k in range(n + 1):
            if k % 2 == 0:                           # ∫(x-c)^k dx = 2h^{k+1}/(k+1)
                coef = fX[k] if k == n else fc[k]
                total = total + coef * (I(w) / I(float(k + 1)))
            w = w * h
    return total

for m, n in [(4, 8), (8, 12), (16, 16)]:
    J = integrate01(m, n)
    print(m, n, repr(J.lo), repr(J.hi), f"{J.hi - J.lo:.1e}")

出力は次のとおりです。

4 8 0.7468241323952737 0.7468241331857042 7.9e-10
8 12 0.74682413281242 0.7468241328124348 1.5e-14
16 16 0.7468241328124131 0.7468241328124405 2.7e-14

真値 $\frac{\sqrt\pi}2\mathrm{erf}(1)=0.746824132812427\ldots$ はどの区間にも入っています。$(m,n)=(8,12)$ で幅は $1.5\times10^{-14}$ になり、それ以上細かくしても幅は縮みません。打ち切りの剰余より、丸めを外へ広げた分の累積が支配的になるためです。計算に使ったのは四則と比較だけなので、出力の区間は第 1 章の 補題 4.7［浮動小数と区間演算］ と本章の定理から、$\int_0^1e^{-x^2}dx$ を含むことが証明されています。

6. 演習

演習 6.1標準

$f(x)=1/(2-x)$ について、(1) $f$ が正則でいられる最大の Bernstein 楕円の $\rho$ を求めよ。(2) $F(w)=f(J(w))$ の Laurent 展開を部分分数で求め、Chebyshev 係数が $a_0=1/\sqrt3$, $a_k=\frac2{\sqrt3}(2-\sqrt3)^k$（$k\ge1$）であることを示せ。

解答

(1) 極は $x=2$ です。$\partial E_\rho$ の長半径が $2$ になる $\rho$ は $\tfrac12(\rho+\rho^{-1})=2$ から $\rho=2+\sqrt3$ です。

(2) $\alpha=2-\sqrt3$, $\beta=2+\sqrt3$（$\alpha\beta=1$, $\beta-\alpha=2\sqrt3$）とおくと、$F(w)=\dfrac{1}{2-\frac12(w+w^{-1})}=\dfrac{-2w}{(w-\alpha)(w-\beta)}=\dfrac{\alpha/\sqrt3}{w-\alpha}-\dfrac{\beta/\sqrt3}{w-\beta}$ です。円環 $\alpha\lt|w|\lt\beta$ では

$$\frac{\alpha/\sqrt3}{w-\alpha}=\frac1{\sqrt3}\sum_{m\ge1}\alpha^mw^{-m},\qquad-\frac{\beta/\sqrt3}{w-\beta}=\frac1{\sqrt3}\sum_{m\ge0}\beta^{-m}w^m$$

なので、$\beta^{-1}=\alpha$ から $c_k=\alpha^{|k|}/\sqrt3$ です。定理 3.1 の証明により $a_0=c_0=1/\sqrt3$, $a_k=2c_k=\frac2{\sqrt3}\alpha^k$ です。検算として $x=1$ で $\sum_ka_k=\frac1{\sqrt3}+\frac2{\sqrt3}\cdot\frac{\alpha}{1-\alpha}=\frac1{\sqrt3}+\frac{\sqrt3-1}{\sqrt3}=1=f(1)$ です（$\frac{\alpha}{1-\alpha}=\frac{\sqrt3-1}2$）。減衰率 $\alpha=1/\rho$ は、定理 3.1 で $r\uparrow\rho$ としたときの $r^{-k}$ と一致します。

演習 6.2易

$f(c)\ne0$ のとき、$g=1/f$ の Taylor 係数が $(g)_0=1/(f)_0$, $(g)_k=-\frac1{(f)_0}\sum_{j=1}^k(f)_j(g)_{k-j}$（$k\ge1$）で与えられることを示せ。

解答

$fg=1$ なので、補題 5.1 の第 1 式から $\sum_{j=0}^k(f)_j(g)_{k-j}$ は $k=0$ で $1$、$k\ge1$ で $0$ です。$k=0$ から $(g)_0=1/(f)_0$、$k\ge1$ では $j=0$ の項 $(f)_0(g)_k$ を分離して解けば主張の式です。区間で評価するときは、$(f)_0$ の区間が $0$ を含まないことを確かめる必要があります。

演習 6.3標準

定理 5.3 を $n=2$, $c=h=\tfrac12$ で $f(x)=e^{-x^2}$ に使い、$\int_0^1e^{-x^2}dx$ を含む区間を手計算で求めよ。$f''$ の $[0,1]$ での値域は、$f'''$ の符号から求めよ。

解答

$G_0\ni f(\tfrac12)=e^{-1/4}\approx0.7788$ です。$f''(x)=(4x^2-2)e^{-x^2}$、$f'''(x)=4x(3-2x^2)e^{-x^2}\ge0$（$0\le x\le1$）なので $f''$ は増加で、値域は $[f''(0),f''(1)]=[-2,\,2e^{-1}]\approx[-2,0.7358]$、$f''/2!$ の値域は $[-1,0.3679]$ です。$\frac{2h^3}3=\frac1{12}$ なので

$$\int_0^1e^{-x^2}dx\in0.7788+[-1,0.3679]\cdot\tfrac1{12}\approx[0.6955,\ 0.8095].$$

真値 $0.7468$ を含みます。剰余の区間が非対称なのは $f''$ が $[0,1]$ でほぼ負であることを反映しており、中点則が過大評価になることも読み取れます。

参考文献

• L. N. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2013 — 第 3 章（Chebyshev 多項式と級数）、第 8 章（解析関数の Chebyshev 係数と Bernstein 楕円）、第 19 章（Gauss 求積と Clenshaw–Curtis 求積の誤差）。
• P. J. Davis, P. Rabinowitz, Methods of Numerical Integration, 2nd ed., Academic Press, 1984 — 第 2 章（Gauss 型公式）、第 4 章（誤差解析）。
• W. Tucker, Validated Numerics: A Short Introduction to Rigorous Computations, Princeton University Press, 2011 — 自動微分と Taylor 係数、検証つき積分。
• R. E. Moore, R. B. Kearfott, M. J. Cloud, Introduction to Interval Analysis, SIAM, 2009 — Taylor 係数の漸化式と区間積分。
• L. B. Rall, Automatic Differentiation: Techniques and Applications, Lecture Notes in Computer Science 120, Springer, 1981.