§E20.20数値積分

最終更新

閉区間上の関数の Riemann 積分は、区間の分割を細かくしていく極限として定まる。この定義は積分の値を有限回の計算で与えるものではないので、積分の値を実際に求めるには、有限個の点での関数値から積分を近似する方法が必要になる。有限個の節点での関数値に重みを掛けて加えた和を求積公式という。節点での値から補間多項式を作り、その積分を近似値とすれば求積公式が得られる。端点の2点による補間から台形則が、中点を加えた3点による補間から Simpson 則が得られ、区間を小区間に分けてこれらを並べたものが複合台形則と複合 Simpson 則である。

求積公式の良さは、何次までの多項式の積分を正確に与えるか（代数的精度）と、滑らかな関数に対する誤差の大きさとで測られる。節点がnn個の求積公式の代数的精度は2n−12n-1を超えることができず、節点を Legendre 多項式の零点にとる Gauss–Legendre 公式はこの上限に達する。たとえば区間[−1,1][-1,1]上で、2点の値の和f(−1/3)+f(1/3)f(-1/\sqrt3)+f(1/\sqrt3)は3次以下のすべての多項式の積分を正確に与える。より一般に、有界開区間上の重み関数に関する直交多項式の零点を節点とする Gauss 求積公式の重みはすべて正であり、このことから、連続関数の重み付き積分に対する誤差は節点の個数を増やすと00に収束する。

誤差の形がわかれば、それを近似の改良に用いることもできる。関数が十分に滑らかならば、複合台形則の誤差は刻み幅の偶数乗による展開をもち、これに Richardson 外挿を繰り返して Romberg 積分が得られる。また、Chebyshev 節点での補間を積分する求積公式は、Bernstein 楕円を含む開集合上で正則な関数に対して、節点の個数について幾何的に減少する誤差の上界をもつ。本記事では、補間にもとづく求積公式の構成と、その誤差の基本的な性質について解説する。

1 補間型求積公式

定義 1.1.a<ba<bを実数とする。

  1. n∈N≥0n\in\Nとし、[a,b][a,b]の相異なる点x0,…,xnx_0,\dots,x_nと実数v0,…,vnv_0,\dots,v_nをとる。[a,b][a,b]上の実数値関数ffにQ(f):=∑i=0nvif(xi)Q(f):=\sum_{i=0}^nv_if(x_i)を対応させる写像QQを、節点x0,…,xnx_0,\dots,x_n、重みv0,…,vnv_0,\dots,v_nの 求積公式 (quadrature rule) といい、viv_iをxix_iにおけるQQの 重み (weight) という。
  2. 節点x0,…,xnx_0,\dots,x_n、重みv0,…,vnv_0,\dots,v_nの求積公式QQが、x0,…,xnx_0,\dots,x_nの Lagrange 基底ℓ0,…,ℓn\ell_0,\dots,\ell_nについてvi=∫abℓi(x) dxv_i=\int_a^b\ell_i(x)\,dx（0≤i≤n0\le i\le n）を満たすとき、QQを 補間型求積公式 (interpolatory quadrature rule) という。n∈N≥1n\in\NNとし、xi=a+i(b−a)/nx_i=a+i(b-a)/n（0≤i≤n0\le i\le n）を節点とする補間型求積公式を、[a,b][a,b]のn+1n+1点の Newton–Cotes 公式 (Newton–Cotes formula) という。
  3. d∈N≥0d\in\Nとし、QQを[a,b][a,b]上の求積公式とする。QQが任意のq∈Pdq\in\mathcal P_dに対してQ(q)=∫abq(x) dxQ(q)=\int_a^bq(x)\,dxを満たすとき、QQの代数的精度はdd以上であるという。QQの代数的精度がdd以上であり、d+1d+1以上でないとき、QQの 代数的精度 (degree of exactness) はddであるという。
  4. c:=(a+b)/2c:=(a+b)/2と置く。[a,b][a,b]上の実数値関数ffに対して T[a,b](f):=b−a2(f(a)+f(b)),S[a,b](f):=b−a6(f(a)+4f(c)+f(b))T_{[a,b]}(f):=\frac{b-a}2\bigl(f(a)+f(b)\bigr),\qquad S_{[a,b]}(f):=\frac{b-a}6\bigl(f(a)+4f(c)+f(b)\bigr) と置き、T[a,b]T_{[a,b]}を[a,b][a,b]の 台形則 (trapezoidal rule)、S[a,b]S_{[a,b]}を[a,b][a,b]の Simpson 則 (Simpson's rule) という。
  5. N∈N≥1N\in\NNとし、h:=(b−a)/Nh:=(b-a)/N、xj:=a+jhx_j:=a+jh（0≤j≤N0\le j\le N）と置く。[a,b][a,b]上の実数値関数ffに対して Th(f):=∑j=0N−1T[xj,xj+1](f)=h(f(x0)2+∑j=1N−1f(xj)+f(xN)2)T_h(f):=\sum_{j=0}^{N-1}T_{[x_j,x_{j+1}]}(f)=h\Bigl(\frac{f(x_0)}2+\sum_{j=1}^{N-1}f(x_j)+\frac{f(x_N)}2\Bigr) を 複合台形則 (composite trapezoidal rule) という。NNが偶数のとき、 Sh(f):=∑k=0N/2−1S[x2k,x2k+2](f)=h3(f(x0)+4∑k=1N/2f(x2k−1)+2∑k=1N/2−1f(x2k)+f(xN))S_h(f):=\sum_{k=0}^{N/2-1}S_{[x_{2k},x_{2k+2}]}(f)=\frac h3\Bigl(f(x_0)+4\sum_{k=1}^{N/2}f(x_{2k-1})+2\sum_{k=1}^{N/2-1}f(x_{2k})+f(x_N)\Bigr) を 複合 Simpson 則 (composite Simpson's rule) という。

補題 1.2.a<ba<bを実数とし、u ⁣:[a,b]→Ru\colon[a,b]\to\Rを連続関数とする。任意のx∈[a,b]x\in[a,b]に対してu(x)≥0u(x)\ge0であり、u(x0)>0u(x_0)>0を満たすx0∈[a,b]x_0\in[a,b]が存在するならば、∫abu(x) dx>0\int_a^bu(x)\,dx>0である。

証明. 定数関数11は(a,b)(a,b)上の重み関数であり、連続関数uuの広義積分∫abu(x)⋅1 dx\int_a^bu(x)\cdot1\,dxは Riemann 積分∫abu(x) dx\int_a^bu(x)\,dxに等しい。§E20.14 補題 1.3 (2)を重み関数11とh=uh=uに適用して主張を得る。▨

命題 1.3.a<ba<bを実数、n∈N≥0n\in\Nとし、QQを[a,b][a,b]の相異なる節点x0,…,xnx_0,\dots,x_n、重みv0,…,vnv_0,\dots,v_nの求積公式とする。ωn+1(x):=∏i=0n(x−xi)\omega_{n+1}(x):=\prod_{i=0}^n(x-x_i)と置く。

  1. QQの代数的精度がnn以上であることと、QQが補間型であることは同値である。QQが補間型ならば、[a,b][a,b]上の任意の実数値関数ffについて、x0,…,xnx_0,\dots,x_nにおけるffの補間多項式ppはQ(f)=∫abp(x) dxQ(f)=\int_a^bp(x)\,dxを満たす。
  2. Q(ωn+12)=0Q(\omega_{n+1}^2)=0かつ∫abωn+1(x)2 dx>0\int_a^b\omega_{n+1}(x)^2\,dx>0である。特に、QQの代数的精度は2n+22n+2以上でない。
  3. T[a,b]T_{[a,b]}は[a,b][a,b]の22点の Newton–Cotes 公式である。∫ab(x−a)(x−b) dx=−(b−a)3/6\int_a^b(x-a)(x-b)\,dx=-(b-a)^3/6であり、T[a,b]T_{[a,b]}は多項式(x−a)(x−b)(x-a)(x-b)に値00を与える。T[a,b]T_{[a,b]}の代数的精度は11である。
  4. c:=(a+b)/2c:=(a+b)/2と置く。S[a,b]S_{[a,b]}は[a,b][a,b]の33点の Newton–Cotes 公式である。∫ab(x−a)(x−c)2(x−b) dx=−(b−a)5/120\int_a^b(x-a)(x-c)^2(x-b)\,dx=-(b-a)^5/120であり、S[a,b]S_{[a,b]}は多項式(x−a)(x−c)2(x−b)(x-a)(x-c)^2(x-b)に値00を与える。S[a,b]S_{[a,b]}の代数的精度は33である。

証明.x0,…,xnx_0,\dots,x_nの Lagrange 基底をℓ0,…,ℓn\ell_0,\dots,\ell_nとする。§E20.12 定理 1.3により、Pn\mathcal P_nの元qqはq=∑i=0nq(xi)ℓiq=\sum_{i=0}^nq(x_i)\ell_iを満たし、ℓi(xj)\ell_i(x_j)はj=ij=iのとき11、j≠ij\ne iのとき00である。

(1)を示す。QQの代数的精度がnn以上ならば、ℓi∈Pn\ell_i\in\mathcal P_nであるから∫abℓi(x) dx=Q(ℓi)=∑jvjℓi(xj)=vi\int_a^b\ell_i(x)\,dx=Q(\ell_i)=\sum_jv_j\ell_i(x_j)=v_iであり、QQは補間型である。QQが補間型であるとし、ffを[a,b][a,b]上の実数値関数、p=∑if(xi)ℓip=\sum_if(x_i)\ell_iをその補間多項式とすると、∫abp(x) dx=∑if(xi)vi=Q(f)\int_a^bp(x)\,dx=\sum_if(x_i)v_i=Q(f)である。q∈Pnq\in\mathcal P_nは自身の補間多項式であるから、Q(q)=∫abq(x) dxQ(q)=\int_a^bq(x)\,dxであり、QQの代数的精度はnn以上である。

(2)を示す。ωn+12\omega_{n+1}^2は各節点で00になるからQ(ωn+12)=0Q(\omega_{n+1}^2)=0である。ωn+12\omega_{n+1}^2は[a,b][a,b]上で連続かつ非負であり、節点でない点で正であるから、補題 1.2により∫abωn+1(x)2 dx>0\int_a^b\omega_{n+1}(x)^2\,dx>0である。ωn+12∈P2n+2\omega_{n+1}^2\in\mathcal P_{2n+2}であるから、QQの代数的精度は2n+22n+2以上でない。

(3)を示す。節点a,ba,bの Lagrange 基底は(b−x)/(b−a)(b-x)/(b-a)と(x−a)/(b−a)(x-a)/(b-a)であり、それぞれの[a,b][a,b]上の積分は(b−a)/2(b-a)/2である。したがってT[a,b]T_{[a,b]}は22点の Newton–Cotes 公式であり、(1)によりその代数的精度は11以上である。t:=x−at:=x-aと置換して

∫ab(x−a)(x−b) dx=∫0b−at(t−(b−a)) dt=(b−a)33−(b−a)32=−(b−a)36\int_a^b(x-a)(x-b)\,dx=\int_0^{b-a}t\bigl(t-(b-a)\bigr)\,dt=\frac{(b-a)^3}3-\frac{(b-a)^3}2=-\frac{(b-a)^3}6

を得る。(x−a)(x−b)(x-a)(x-b)はa,ba,bで00になるのでT[a,b]T_{[a,b]}はこれに値00を与え、(x−a)(x−b)∈P2(x-a)(x-b)\in\mathcal P_2であるから、T[a,b]T_{[a,b]}の代数的精度は11である。

(4)を示す。h:=(b−a)/2h:=(b-a)/2とし、t:=(x−c)/ht:=(x-c)/hと置換する。節点a,c,ba,c,bの Lagrange 基底はt(t−1)/2t(t-1)/2、1−t21-t^2、t(t+1)/2t(t+1)/2であり、∫ab dx=h∫−11 dt\int_a^b\,dx=h\int_{-1}^1\,dtから、それぞれの積分はh/3=(b−a)/6h/3=(b-a)/6、4h/3=4(b−a)/64h/3=4(b-a)/6、h/3=(b−a)/6h/3=(b-a)/6である。したがってS[a,b]S_{[a,b]}は33点の Newton–Cotes 公式であり、その代数的精度は22以上である。∫ab(x−c)3 dx=h4∫−11t3 dt=0\int_a^b(x-c)^3\,dx=h^4\int_{-1}^1t^3\,dt=0であり、S[a,b]((x−c)3)=b−a6((−h)3+0+h3)=0S_{[a,b]}\bigl((x-c)^3\bigr)=\frac{b-a}6\bigl((-h)^3+0+h^3\bigr)=0である。P3\mathcal P_3の元はP2\mathcal P_2の元と(x−c)3(x-c)^3の定数倍の和であるから、S[a,b]S_{[a,b]}の代数的精度は33以上である。(x−a)(x−c)2(x−b)=h4(t4−t2)(x-a)(x-c)^2(x-b)=h^4(t^4-t^2)であるから

∫ab(x−a)(x−c)2(x−b) dx=h5(25−23)=−4h515=−(b−a)5120\int_a^b(x-a)(x-c)^2(x-b)\,dx=h^5\Bigl(\frac25-\frac23\Bigr)=-\frac{4h^5}{15}=-\frac{(b-a)^5}{120}

である。この多項式はa,c,ba,c,bで00になるのでS[a,b]S_{[a,b]}はこれに値00を与え、この多項式はP4\mathcal P_4に属するから、S[a,b]S_{[a,b]}の代数的精度は33である。▨

2 台形則と Simpson 則の誤差

補題 2.1.a<ba<b、c<dc<dを実数とし、wwを(a,b)(a,b)上の重み関数とする。K ⁣:[a,b]→RK\colon[a,b]\to\Rは連続であり、すべてのx∈[a,b]x\in[a,b]でK(x)≥0K(x)\ge0であるか、すべてのx∈[a,b]x\in[a,b]でK(x)≤0K(x)\le0であるとする。G ⁣:[c,d]→RG\colon[c,d]\to\Rを連続関数とする。

  1. 連続関数E ⁣:[a,b]→RE\colon[a,b]\to\Rと、各x∈[a,b]x\in[a,b]に対する点ηx∈[c,d]\eta_x\in[c,d]が、すべてのx∈[a,b]x\in[a,b]でE(x)=G(ηx)K(x)E(x)=G(\eta_x)K(x)を満たすならば、 ∫abE(x)w(x) dx=G(ξ)∫abK(x)w(x) dx\int_a^bE(x)w(x)\,dx=G(\xi)\int_a^bK(x)w(x)\,dx を満たすξ∈[c,d]\xi\in[c,d]が存在する。対応x↦ηxx\mapsto\eta_xには条件を課さない。
  2. (1)の仮定に加えて、すべてのx∈[a,b]x\in[a,b]でηx∈(c,d)\eta_x\in(c,d)ならば、(1)のξ\xiを(c,d)(c,d)にとることができる。
  3. [c,d]=[a,b][c,d]=[a,b]ならば、∫abK(x)G(x)w(x) dx=G(ξ)∫abK(x)w(x) dx\int_a^bK(x)G(x)w(x)\,dx=G(\xi)\int_a^bK(x)w(x)\,dxを満たすξ∈[a,b]\xi\in[a,b]が存在する。

w≡1w\equiv1とすると、各積分は[a,b][a,b]上の Riemann 積分であり、上の三つの主張は Riemann 積分についての主張として成り立つ。

証明.[a,b][a,b]上の連続関数ggについて、広義積分∫abg(x)w(x) dx\int_a^bg(x)w(x)\,dxは§E20.14 補題 1.3 (1)により定まり、w≡1w\equiv1のときは Riemann 積分∫abg(x) dx\int_a^bg(x)\,dxに等しい。g≥0g\ge0ならば、ggが恒等的に00のとき∫abgw dx=0\int_a^bgw\,dx=0であり、そうでないとき§E20.14 補題 1.3 (2)により∫abgw dx>0\int_a^bgw\,dx>0である。

(1)と(2)を示す。K,EK,Eを−K,−E-K,-Eに置き換えても仮定と結論は変わらないから、K≥0K\ge0と仮定する。§D1.13 定理 2.1によりG(y−)=m:=min⁡[c,d]GG(y_-)=m:=\min_{[c,d]}G、G(y+)=M:=max⁡[c,d]GG(y_+)=M:=\max_{[c,d]}Gを満たすy±∈[c,d]y_\pm\in[c,d]が存在する。すべてのxxでmK(x)≤G(ηx)K(x)=E(x)≤MK(x)mK(x)\le G(\eta_x)K(x)=E(x)\le MK(x)である。KKが恒等的に00ならばEEも恒等的に00であり、任意のξ∈(c,d)\xi\in(c,d)で等式が成り立つ。KKが恒等的に00でないとき、IK:=∫abKw dx>0I_K:=\int_a^bKw\,dx>0であり、連続な非負関数E−mKE-mKとMK−EMK-Eの積分は00以上であるから、μ:=∫abEw dx/IK\mu:=\int_a^bEw\,dx/I_Kはm≤μ≤Mm\le\mu\le Mを満たす。m<μ<Mm<\mu<Mならば、§D1.12 系 1.2をy−y_-とy+y_+を端点とする閉区間上のGGに適用してG(ξ)=μG(\xi)=\muを満たすξ\xiを得る。G(ξ)≠m,MG(\xi)\ne m,Mからξ≠y±\xi\ne y_\pmであり、ξ\xiはy−y_-とy+y_+の間の開区間に属するから、ξ∈(c,d)\xi\in(c,d)である。μ=m\mu=mならば、∫ab(E−mK)w dx=0\int_a^b(E-mK)w\,dx=0であるからE−mKE-mKは恒等的に00である。K(x0)>0K(x_0)>0を満たすx0x_0をとるとG(ηx0)K(x0)=mK(x0)G(\eta_{x_0})K(x_0)=mK(x_0)であり、ξ:=ηx0\xi:=\eta_{x_0}はG(ξ)=m=μG(\xi)=m=\muを満たす。μ=M\mu=Mならば、MK−EMK-Eに同じ議論を適用して、ξ:=ηx0\xi:=\eta_{x_0}はG(ξ)=M=μG(\xi)=M=\muを満たす。後の二つの場合のξ\xiはηx0\eta_{x_0}であるから、(2)の仮定の下で(c,d)(c,d)に属する。

(3)は、(1)をηx:=x\eta_x:=x、E:=KGE:=KGとして適用したものである。▨

補題 2.2.a<ba<bを実数、N∈N≥1N\in\NNとし、g ⁣:[a,b]→Rg\colon[a,b]\to\Rを連続関数、η1,…,ηN∈(a,b)\eta_1,\dots,\eta_N\in(a,b)とする。このときg(ξ)=1N∑j=1Ng(ηj)g(\xi)=\frac1N\sum_{j=1}^Ng(\eta_j)を満たすξ∈(a,b)\xi\in(a,b)が存在する。

証明.g(ηp)=min⁡jg(ηj)g(\eta_p)=\min_jg(\eta_j)、g(ηq)=max⁡jg(ηj)g(\eta_q)=\max_jg(\eta_j)を満たすp,qp,qをとる。平均μ:=1N∑jg(ηj)\mu:=\frac1N\sum_jg(\eta_j)はg(ηp)≤μ≤g(ηq)g(\eta_p)\le\mu\le g(\eta_q)を満たす。ηp=ηq\eta_p=\eta_qならばμ=g(ηp)\mu=g(\eta_p)であり、ξ:=ηp\xi:=\eta_pとする。ηp≠ηq\eta_p\ne\eta_qならば、§D1.12 系 1.2をηp\eta_pとηq\eta_qを端点とする閉区間上のggに適用してg(ξ)=μg(\xi)=\muを満たすξ\xiを得る。この閉区間は(a,b)(a,b)に含まれる。▨

定理 2.3.a<ba<bを実数とし、f ⁣:[a,b]→Rf\colon[a,b]\to\RをC2C^2級の関数とする。

  1. ∫abf(x) dx−T[a,b](f)=−(b−a)312f′′(ξ)\int_a^bf(x)\,dx-T_{[a,b]}(f)=-\dfrac{(b-a)^3}{12}f''(\xi)を満たすξ∈(a,b)\xi\in(a,b)が存在する。
  2. N∈N≥1N\in\NN、h:=(b−a)/Nh:=(b-a)/Nならば、∫abf(x) dx−Th(f)=−(b−a)h212f′′(ξ)\int_a^bf(x)\,dx-T_h(f)=-\dfrac{(b-a)h^2}{12}f''(\xi)を満たすξ∈(a,b)\xi\in(a,b)が存在する。

証明.(1)を示す。ppをa,ba,bにおけるffの補間多項式とする。命題 1.3 (3)と命題 1.3 (1)によりT[a,b](f)=∫abp(x) dxT_{[a,b]}(f)=\int_a^bp(x)\,dxである。§E20.12 定理 2.2 (1)をn=1n=1として適用すると、各x∈[a,b]x\in[a,b]に対してf(x)−p(x)=f′′(ηx)2(x−a)(x−b)f(x)-p(x)=\frac{f''(\eta_x)}2(x-a)(x-b)を満たすηx∈(a,b)\eta_x\in(a,b)が存在する。E:=f−pE:=f-pは連続であり、K(x):=(x−a)(x−b)/2K(x):=(x-a)(x-b)/2は[a,b][a,b]上で00以下である。補題 2.1 (2)をw≡1w\equiv1、[c,d]=[a,b][c,d]=[a,b]、G=f′′G=f''として適用し、命題 1.3 (3)の積分の値を用いると、

∫abf(x) dx−T[a,b](f)=∫abE(x) dx=f′′(ξ)∫ab(x−a)(x−b)2 dx=−(b−a)312f′′(ξ)\int_a^bf(x)\,dx-T_{[a,b]}(f)=\int_a^bE(x)\,dx=f''(\xi)\int_a^b\frac{(x-a)(x-b)}2\,dx=-\frac{(b-a)^3}{12}f''(\xi)

を満たすξ∈(a,b)\xi\in(a,b)を得る。

(2)を示す。xj:=a+jhx_j:=a+jhと置く。(1)を各[xj,xj+1][x_j,x_{j+1}]上のffに適用して、∫xjxj+1f(x) dx−T[xj,xj+1](f)=−h312f′′(ξj)\int_{x_j}^{x_{j+1}}f(x)\,dx-T_{[x_j,x_{j+1}]}(f)=-\frac{h^3}{12}f''(\xi_j)を満たすξj∈(xj,xj+1)⊆(a,b)\xi_j\in(x_j,x_{j+1})\subseteq(a,b)をとる。§D1.17 定理 3.6とThT_hの定義によりjjについて和をとり、Nh=b−aNh=b-aを用いると

∫abf(x) dx−Th(f)=−Nh312⋅1N∑j=0N−1f′′(ξj)=−(b−a)h212⋅1N∑j=0N−1f′′(ξj)\int_a^bf(x)\,dx-T_h(f)=-\frac{Nh^3}{12}\cdot\frac1N\sum_{j=0}^{N-1}f''(\xi_j)=-\frac{(b-a)h^2}{12}\cdot\frac1N\sum_{j=0}^{N-1}f''(\xi_j)

である。補題 2.2をg=f′′g=f''に適用して主張を得る。▨

定理 2.4.a<ba<bを実数とし、f ⁣:[a,b]→Rf\colon[a,b]\to\RをC4C^4級の関数とする。

  1. h:=(b−a)/2h:=(b-a)/2と置くと、∫abf(x) dx−S[a,b](f)=−h590f(4)(ξ)\int_a^bf(x)\,dx-S_{[a,b]}(f)=-\dfrac{h^5}{90}f^{(4)}(\xi)を満たすξ∈(a,b)\xi\in(a,b)が存在する。
  2. NNが正の偶数であり、h:=(b−a)/Nh:=(b-a)/Nならば、∫abf(x) dx−Sh(f)=−(b−a)h4180f(4)(ξ)\int_a^bf(x)\,dx-S_h(f)=-\dfrac{(b-a)h^4}{180}f^{(4)}(\xi)を満たすξ∈(a,b)\xi\in(a,b)が存在する。

証明.(1)を示す。c:=(a+b)/2c:=(a+b)/2、Ω3(x):=(x−a)(x−c)(x−b)\Omega_3(x):=(x-a)(x-c)(x-b)、Ω(x):=(x−a)(x−c)2(x−b)\Omega(x):=(x-a)(x-c)^2(x-b)と置く。ppをa,c,ba,c,bにおけるffの補間多項式とし、

p3:=p+p′(c)−f′(c)h2 Ω3p_3:=p+\frac{p'(c)-f'(c)}{h^2}\,\Omega_3

と置く。Ω3′(c)=(c−a)(c−b)=−h2\Omega_3'(c)=(c-a)(c-b)=-h^2であるから、p3∈P3p_3\in\mathcal P_3はa,c,ba,c,bでffと同じ値をとり、p3′(c)=f′(c)p_3'(c)=f'(c)を満たす。S[a,b](f)=S[a,b](p3)S_{[a,b]}(f)=S_{[a,b]}(p_3)であり、命題 1.3 (4)によりS[a,b](p3)=∫abp3(x) dxS_{[a,b]}(p_3)=\int_a^bp_3(x)\,dxである。

x∈[a,b]x\in[a,b]とする。x∈{a,c,b}x\in\{a,c,b\}ならばf(x)−p3(x)=0=Ω(x)f(x)-p_3(x)=0=\Omega(x)である。x∉{a,c,b}x\notin\{a,c,b\}ならば、κ:=(f(x)−p3(x))/Ω(x)\kappa:=(f(x)-p_3(x))/\Omega(x)、g(t):=f(t)−p3(t)−κΩ(t)g(t):=f(t)-p_3(t)-\kappa\Omega(t)（t∈[a,b]t\in[a,b]）と置く。ggはC4C^4級であり、Ω\Omegaは因子(t−c)2(t-c)^2をもつからΩ(c)=Ω′(c)=0\Omega(c)=\Omega'(c)=0であって、g(a)=g(b)=g(c)=g′(c)=0g(a)=g(b)=g(c)=g'(c)=0、g(x)=0g(x)=0である。相異なる四点a,c,b,xa,c,b,xに重複度1,2,1,11,2,1,1を与えると、重複度は11以上44以下でその和は55であり、四点の最小点はaa、最大点はbbである。§E20.12 補題 2.1をm=4m=4として適用して、g(4)(ηx)=0g^{(4)}(\eta_x)=0を満たすηx∈(a,b)\eta_x\in(a,b)を得る。p3(4)=0p_3^{(4)}=0、Ω(4)=24\Omega^{(4)}=24であるからκ=f(4)(ηx)/24\kappa=f^{(4)}(\eta_x)/24である。したがって、各x∈[a,b]x\in[a,b]に対して

f(x)−p3(x)=f(4)(ηx)24 Ω(x)f(x)-p_3(x)=\frac{f^{(4)}(\eta_x)}{24}\,\Omega(x)

を満たすηx∈(a,b)\eta_x\in(a,b)が存在する。

f−p3f-p_3は連続であり、Ω\Omegaは[a,b][a,b]上で00以下である。補題 2.1 (2)を、w≡1w\equiv1、K=ΩK=\Omega、[a,b][a,b]上の関数G=f(4)/24G=f^{(4)}/24として適用し、命題 1.3 (4)の積分の値と(b−a)5=32h5(b-a)^5=32h^5を用いると、

∫abf(x) dx−S[a,b](f)=∫ab(f(x)−p3(x)) dx=f(4)(ξ)24⋅(−(b−a)5120)=−h590f(4)(ξ)\int_a^bf(x)\,dx-S_{[a,b]}(f)=\int_a^b\bigl(f(x)-p_3(x)\bigr)\,dx=\frac{f^{(4)}(\xi)}{24}\cdot\Bigl(-\frac{(b-a)^5}{120}\Bigr)=-\frac{h^5}{90}f^{(4)}(\xi)

を満たすξ∈(a,b)\xi\in(a,b)を得る。

(2)を示す。xj:=a+jhx_j:=a+jhと置く。0≤k≤N/2−10\le k\le N/2-1について、(1)を長さ2h2hの区間[x2k,x2k+2][x_{2k},x_{2k+2}]上のffに適用して、∫x2kx2k+2f(x) dx−S[x2k,x2k+2](f)=−h590f(4)(ξk)\int_{x_{2k}}^{x_{2k+2}}f(x)\,dx-S_{[x_{2k},x_{2k+2}]}(f)=-\frac{h^5}{90}f^{(4)}(\xi_k)を満たすξk∈(x2k,x2k+2)⊆(a,b)\xi_k\in(x_{2k},x_{2k+2})\subseteq(a,b)をとる。§D1.17 定理 3.6とShS_hの定義によりkkについて和をとり、(N/2)⋅h=(b−a)/2(N/2)\cdot h=(b-a)/2を用いると

∫abf(x) dx−Sh(f)=−(b−a)h4180⋅2N∑k=0N/2−1f(4)(ξk)\int_a^bf(x)\,dx-S_h(f)=-\frac{(b-a)h^4}{180}\cdot\frac2N\sum_{k=0}^{N/2-1}f^{(4)}(\xi_k)

である。補題 2.2をg=f(4)g=f^{(4)}とN/2N/2個の点ξk\xi_kに適用して主張を得る。▨

例 2.5.f(x):=exf(x):=e^x、[a,b]=[0,1][a,b]=[0,1]、N=4N=4、h=1/4h=1/4とする。Th(f)T_h(f)とSh(f)S_h(f)は、ともに同じ55個の値f(0),f(1/4),f(1/2),f(3/4),f(1)f(0),f(1/4),f(1/2),f(3/4),f(1)から計算される。5050桁の十進演算で計算すると、∫01ex dx=e−1=1.718281828459…\int_0^1e^x\,dx=e-1=1.718281828459\ldotsに対して

Th(f)=1.727221904557…,Sh(f)=1.718318841921…T_h(f)=1.727221904557\ldots,\qquad S_h(f)=1.718318841921\ldots

であり、誤差は∫01f−Th(f)=−8.940076098…×10−3\int_0^1f-T_h(f)=-8.940076098\ldots\times10^{-3}、∫01f−Sh(f)=−3.701346270…×10−5\int_0^1f-S_h(f)=-3.701346270\ldots\times10^{-5}である。f′′=f(4)=exf''=f^{(4)}=e^xは(0,1)(0,1)上で11とeeの間の値をとるから、定理 2.3 (2)により∫01f−Th(f)=−eξT/192\int_0^1f-T_h(f)=-e^{\xi_T}/192、定理 2.4 (2)により∫01f−Sh(f)=−eξS/46080\int_0^1f-S_h(f)=-e^{\xi_S}/46080を満たすξT,ξS∈(0,1)\xi_T,\xi_S\in(0,1)が存在し、

−e192<∫01f−Th(f)<−1192,−e46080<∫01f−Sh(f)<−146080-\frac e{192}<\int_0^1f-T_h(f)<-\frac1{192},\qquad-\frac e{46080}<\int_0^1f-S_h(f)<-\frac1{46080}

である。数値で書くと、第一の区間は(−1.4158×10−2,−5.2083×10−3)(-1.4158\times10^{-2},-5.2083\times10^{-3})、第二の区間は(−5.8990×10−5,−2.1701×10−5)(-5.8990\times10^{-5},-2.1701\times10^{-5})であり、計算した誤差はそれぞれの区間に属する。二つの誤差の比は240 eξT−ξS240\,e^{\xi_T-\xi_S}に等しく、240=15/h2240=15/h^2は二つの誤差公式の係数(b−a)h2/12(b-a)h^2/12と(b−a)h4/180(b-a)h^4/180の比である。計算した誤差の比は241.53…241.53\ldotsである。

例 2.6.f(x):=∣x−1/2∣f(x):=\lvert x-1/2\rvert（x∈[0,1]x\in[0,1]）は1/21/2で微分可能でなく、定理 2.4の仮定を満たさない。NNを44で割った余りが22である正の整数とし、h:=1/Nh:=1/N、xj:=jhx_j:=jhと置く。N/2N/2は奇数であるから1/2=xN/21/2=x_{N/2}は添字が奇数の格子点であり、ShS_hを定める区間のうち[1/2−h,1/2+h][1/2-h,1/2+h]の中点である。その他の区間[x2k,x2k+2][x_{2k},x_{2k+2}]ではffは一次式であるから、命題 1.3 (4)により Simpson 則は積分に等しい値を与える。区間[1/2−h,1/2+h][1/2-h,1/2+h]では∫1/2−h1/2+hf(x) dx=h2\int_{1/2-h}^{1/2+h}f(x)\,dx=h^2、S[1/2−h,1/2+h](f)=2h6(h+0+h)=2h23S_{[1/2-h,1/2+h]}(f)=\frac{2h}6(h+0+h)=\frac{2h^2}3である。したがって

∫01f(x) dx−Sh(f)=h23\int_0^1f(x)\,dx-S_h(f)=\frac{h^2}3

である。定数CCがすべての正の偶数NNについて∣∫01f−Sh(f)∣≤Ch4\bigl|\int_0^1f-S_h(f)\bigr|\le Ch^4を満たすならば、N≡2(mod4)N\equiv2\pmod4についてC≥N2/3C\ge N^2/3となり、そのようなCCは存在しない。同じNNについて1/21/2は格子点であり、各[xj,xj+1][x_j,x_{j+1}]でffは一次式であるから、命題 1.3 (3)によりTh(f)=∫01f(x) dx=1/4T_h(f)=\int_0^1f(x)\,dx=1/4である。NNが44の倍数ならば、1/2=xN/21/2=x_{N/2}は区間[x2k,x2k+2][x_{2k},x_{2k+2}]の端点であり、Sh(f)=∫01f(x) dxS_h(f)=\int_0^1f(x)\,dxである。

3 Gauss 求積

定義 3.1.a<ba<bを実数、wwを(a,b)(a,b)上の重み関数、n∈N≥1n\in\NNとし、πn\pi_nをwwに関するnn次の直交多項式とする。§E20.14 定理 3.1により、πn=∏i=1n(x−zi)\pi_n=\prod_{i=1}^n(x-z_i)を満たす実数a<z1<⋯<zn<ba<z_1<\dots<z_n<bがただ一つ存在する。1≤i≤n1\le i\le nに対してLi(x):=∏j≠i(x−zj)/(zi−zj)L_i(x):=\prod_{j\ne i}(x-z_j)/(z_i-z_j)（n=1n=1のときL1:=1L_1:=1）、λi:=∫abLi(x)w(x) dx\lambda_i:=\int_a^bL_i(x)w(x)\,dx（§E20.14 補題 1.3 (1)により有限）と置く。節点z1,…,znz_1,\dots,z_n、重みλ1,…,λn\lambda_1,\dots,\lambda_nの求積公式

Gnw(f):=∑i=1nλif(zi)G_n^w(f):=\sum_{i=1}^n\lambda_if(z_i)

を、wwに関するnn点の Gauss 求積公式 (Gaussian quadrature rule) という。

定理 3.2.a<ba<bを実数、wwを(a,b)(a,b)上の重み関数、μ0:=∫abw(x) dx\mu_0:=\int_a^bw(x)\,dx、n∈N≥1n\in\NNとし、πn\pi_n、ziz_i、LiL_i、λi\lambda_i、GnwG_n^wを定義 3.1のとおりとする。

  1. 任意のq∈P2n−1q\in\mathcal P_{2n-1}に対してGnw(q)=∫abq(x)w(x) dxG_n^w(q)=\int_a^bq(x)w(x)\,dxである。またGnw(πn2)=0<∫abπn(x)2w(x) dxG_n^w(\pi_n^2)=0<\int_a^b\pi_n(x)^2w(x)\,dxである。特にw≡1w\equiv1のとき、GnwG_n^wの代数的精度は2n−12n-1である。
  2. 1≤i≤n1\le i\le nに対してλi=∫abLi(x)2w(x) dx>0\lambda_i=\int_a^bL_i(x)^2w(x)\,dx>0であり、∑i=1nλi=μ0\sum_{i=1}^n\lambda_i=\mu_0である。
  3. f ⁣:[a,b]→Rf\colon[a,b]\to\RがC2nC^{2n}級ならば、 ∫abf(x)w(x) dx−Gnw(f)=f(2n)(ξ)(2n)!∫abπn(x)2w(x) dx\int_a^bf(x)w(x)\,dx-G_n^w(f)=\frac{f^{(2n)}(\xi)}{(2n)!}\int_a^b\pi_n(x)^2w(x)\,dx を満たすξ∈(a,b)\xi\in(a,b)が存在する。
  4. f∈C([a,b])f\in C([a,b])とし、∥⋅∥∞\lVert\cdot\rVert_\inftyを[a,b][a,b]上の最大値ノルム、E2n−1(f)E_{2n-1}(f)をffのP2n−1\mathcal P_{2n-1}による最良一様近似誤差とする。任意のp∈P2n−1p\in\mathcal P_{2n-1}に対して ∣∫abf(x)w(x) dx−Gnw(f)∣≤2μ0∥f−p∥∞\Bigl|\int_a^bf(x)w(x)\,dx-G_n^w(f)\Bigr|\le2\mu_0\lVert f-p\rVert_\infty であり、左辺は2μ0E2n−1(f)2\mu_0E_{2n-1}(f)以下である。各n∈N≥1n\in\NNに対するGnwG_n^wについて、左辺はn→∞n\to\inftyのとき00に収束する。

証明.§E20.12 定理 1.3を相異なるnn点z1,…,znz_1,\dots,z_nに適用すると、r∈Pn−1r\in\mathcal P_{n-1}はr=∑ir(zi)Lir=\sum_ir(z_i)L_iを満たし、Li(zj)L_i(z_j)はj=ij=iのとき11、j≠ij\ne iのとき00である。

(1)を示す。q∈P2n−1q\in\mathcal P_{2n-1}をモニック多項式πn\pi_nで割り、q=sπn+rq=s\pi_n+r、s,r∈Pn−1s,r\in\mathcal P_{n-1}と書く。πn(zi)=0\pi_n(z_i)=0であるからGnw(q)=∑iλir(zi)=∫ab∑ir(zi)Li(x)w(x) dx=∫abr(x)w(x) dxG_n^w(q)=\sum_i\lambda_ir(z_i)=\int_a^b\sum_ir(z_i)L_i(x)w(x)\,dx=\int_a^br(x)w(x)\,dxである。直交多項式の定義により∫absπnw dx=⟨πn,s⟩w=0\int_a^bs\pi_nw\,dx=\langle\pi_n,s\rangle_w=0であるから、∫abqw dx=∫abrw dx=Gnw(q)\int_a^bqw\,dx=\int_a^brw\,dx=G_n^w(q)である。πn2\pi_n^2は各ziz_iで00になるからGnw(πn2)=0G_n^w(\pi_n^2)=0である。πn2\pi_n^2は[a,b][a,b]上で連続かつ非負であり、z1,…,znz_1,\dots,z_n以外の点で正であるから、§E20.14 補題 1.3 (2)により∫abπn2w dx>0\int_a^b\pi_n^2w\,dx>0である。w≡1w\equiv1のとき∫abqw dx\int_a^bqw\,dxは Riemann 積分∫abq\int_a^bqであり、πn2∈P2n\pi_n^2\in\mathcal P_{2n}であるから、GnwG_n^wの代数的精度は2n−12n-1である。

(2)を示す。Li2∈P2n−2⊆P2n−1L_i^2\in\mathcal P_{2n-2}\subseteq\mathcal P_{2n-1}であるから、(1)により∫abLi2w dx=Gnw(Li2)=∑jλjLi(zj)2=λi\int_a^bL_i^2w\,dx=G_n^w(L_i^2)=\sum_j\lambda_jL_i(z_j)^2=\lambda_iである。Li2L_i^2は連続かつ非負であり、Li(zi)2=1>0L_i(z_i)^2=1>0であるから、§E20.14 補題 1.3 (2)によりλi>0\lambda_i>0である。定数関数1∈P2n−11\in\mathcal P_{2n-1}に(1)を適用して∑iλi=Gnw(1)=μ0\sum_i\lambda_i=G_n^w(1)=\mu_0を得る。

(3)を示す。HHをz1,…,znz_1,\dots,z_nにおけるffの Hermite 補間多項式とする（§E20.12 定理 7.2）。H(zi)=f(zi)H(z_i)=f(z_i)であるからGnw(H)=Gnw(f)G_n^w(H)=G_n^w(f)であり、H∈P2n−1H\in\mathcal P_{2n-1}であるから(1)によりGnw(H)=∫abHw dxG_n^w(H)=\int_a^bHw\,dxである。したがって∫abfw dx−Gnw(f)=∫ab(f−H)w dx\int_a^bfw\,dx-G_n^w(f)=\int_a^b(f-H)w\,dxである。§E20.12 定理 7.3により、各x∈[a,b]x\in[a,b]に対してf(x)−H(x)=f(2n)(ηx)(2n)!πn(x)2f(x)-H(x)=\frac{f^{(2n)}(\eta_x)}{(2n)!}\pi_n(x)^2を満たすηx∈(a,b)\eta_x\in(a,b)が存在する。f−Hf-Hは連続であり、πn2\pi_n^2は非負である。補題 2.1 (2)を[c,d]=[a,b][c,d]=[a,b]、K=πn2K=\pi_n^2、G=f(2n)/(2n)!G=f^{(2n)}/(2n)!、E=f−HE=f-Hとして適用して主張を得る。

(4)を示す。p∈P2n−1p\in\mathcal P_{2n-1}とする。(1)により∫abpw dx=Gnw(p)\int_a^bpw\,dx=G_n^w(p)であるから、

∫abfw dx−Gnw(f)=∫ab(f−p)w dx−∑i=1nλi(f(zi)−p(zi))\int_a^bfw\,dx-G_n^w(f)=\int_a^b(f-p)w\,dx-\sum_{i=1}^n\lambda_i\bigl(f(z_i)-p(z_i)\bigr)

である。§E20.14 補題 1.3 (1)により第一項の絶対値はμ0∥f−p∥∞\mu_0\lVert f-p\rVert_\infty以下である。第二項の絶対値は∑i∣λi∣ ∥f−p∥∞\sum_i\lvert\lambda_i\rvert\,\lVert f-p\rVert_\infty以下であり、(2)によりすべてのλi\lambda_iは正であるから、∑i∣λi∣=∑iλi=μ0\sum_i\lvert\lambda_i\rvert=\sum_i\lambda_i=\mu_0である。したがって左辺の絶対値は2μ0∥f−p∥∞2\mu_0\lVert f-p\rVert_\infty以下であり、p∈P2n−1p\in\mathcal P_{2n-1}について下限をとると2μ0E2n−1(f)2\mu_0E_{2n-1}(f)以下である。§E20.15 命題 3.1 (1)によりm→∞m\to\inftyのときEm(f)→0E_m(f)\to0であり、2n−1→∞2n-1\to\inftyであるからE2n−1(f)→0E_{2n-1}(f)\to0である。▨

4 Gauss–Legendre 公式

系 4.1.n∈N≥1n\in\NNとし、PnP_nをnn次の Legendre 多項式、cnc_nをPnP_nのxnx^nの係数とする。PnP_nの実数の零点は相異なるnn個の点x1<⋯<xnx_1<\dots<x_nであり、すべて(−1,1)(-1,1)に属し、Pn=cn∏i=1n(x−xi)P_n=c_n\prod_{i=1}^n(x-x_i)である。

証明.(−1,1)(-1,1)上の重み関数11に関するnn次の直交多項式をπn\pi_nとする。§E20.14 命題 4.1によりcn=(2n)!/(2n(n!)2)≠0c_n=(2n)!/(2^n(n!)^2)\ne0かつPn=cnπnP_n=c_n\pi_nである。§E20.14 定理 3.1によりπn=∏i=1n(x−xi)\pi_n=\prod_{i=1}^n(x-x_i)を満たす実数−1<x1<⋯<xn<1-1<x_1<\dots<x_n<1が存在する。したがってPn=cn∏i(x−xi)P_n=c_n\prod_i(x-x_i)であり、cn≠0c_n\ne0であるからPnP_nの実数の零点はx1,…,xnx_1,\dots,x_nである。▨

定義 4.2.n∈N≥1n\in\NNとし、x1<⋯<xnx_1<\dots<x_nをnn次の Legendre 多項式PnP_nの零点（系 4.1）とする。1≤i≤n1\le i\le nに対してLi(x):=∏j≠i(x−xj)/(xi−xj)L_i(x):=\prod_{j\ne i}(x-x_j)/(x_i-x_j)（n=1n=1のときL1:=1L_1:=1）、wi:=∫−11Li(x) dxw_i:=\int_{-1}^1L_i(x)\,dxと置く。

  1. [−1,1][-1,1]上の実数値関数ffに対してGn(f):=∑i=1nwif(xi)G_n(f):=\sum_{i=1}^nw_if(x_i)と置き、GnG_nをnn点の Gauss–Legendre 公式 (Gauss–Legendre quadrature rule) という。
  2. a<ba<bを実数とし、ϕ(t):=a+b2+b−a2t\phi(t):=\frac{a+b}2+\frac{b-a}2tと置く。[a,b][a,b]上の実数値関数ffに対してGn[a,b](f):=b−a2∑i=1nwif(ϕ(xi))G_n^{[a,b]}(f):=\frac{b-a}2\sum_{i=1}^nw_if(\phi(x_i))と置き、Gn[a,b]G_n^{[a,b]}を[a,b][a,b]のnn点の Gauss–Legendre 公式という。

系 4.3.n∈N≥1n\in\NNとし、xix_i、wiw_i、GnG_n、ϕ\phi、Gn[a,b]G_n^{[a,b]}を定義 4.2のとおりとし、πn(x):=∏i=1n(x−xi)\pi_n(x):=\prod_{i=1}^n(x-x_i)と置く。

  1. GnG_nは(−1,1)(-1,1)上の重み関数11に関するnn点の Gauss 求積公式である。したがって、GnG_nの代数的精度は2n−12n-1であり、すべてのiiでwi>0w_i>0、∑i=1nwi=2\sum_{i=1}^nw_i=2である。f ⁣:[−1,1]→Rf\colon[-1,1]\to\RがC2nC^{2n}級ならば、∫−11f(x) dx−Gn(f)=f(2n)(ξ)(2n)!∫−11πn(x)2 dx\int_{-1}^1f(x)\,dx-G_n(f)=\dfrac{f^{(2n)}(\xi)}{(2n)!}\int_{-1}^1\pi_n(x)^2\,dxを満たすξ∈(−1,1)\xi\in(-1,1)が存在する。
  2. a<ba<bを実数とする。Gn[a,b]G_n^{[a,b]}の代数的精度は2n−12n-1である。f ⁣:[a,b]→Rf\colon[a,b]\to\RがC2nC^{2n}級ならば、 ∫abf(x) dx−Gn[a,b](f)=(b−a2)2n+1f(2n)(ξ)(2n)!∫−11πn(x)2 dx\int_a^bf(x)\,dx-G_n^{[a,b]}(f)=\Bigl(\frac{b-a}2\Bigr)^{2n+1}\frac{f^{(2n)}(\xi)}{(2n)!}\int_{-1}^1\pi_n(x)^2\,dx を満たすξ∈(a,b)\xi\in(a,b)が存在する。

証明.(1)を示す。定数関数11は(−1,1)(-1,1)上の重み関数であり、μ0=2\mu_0=2である。§E20.14 命題 4.1と系 4.1により、重み関数11に関するnn次の直交多項式はPn/cn=πnP_n/c_n=\pi_nであるから、GnG_nはその Gauss 求積公式であり、残りの主張は定理 3.2 (1)、定理 3.2 (2)、定理 3.2 (3)を重み関数11に適用したものである。

(2)を示す。[a,b][a,b]上の連続関数ffについて、置換積分x=ϕ(t)x=\phi(t)により∫abf(x) dx=b−a2∫−11f(ϕ(t)) dt\int_a^bf(x)\,dx=\frac{b-a}2\int_{-1}^1f(\phi(t))\,dtであり、定義によりGn[a,b](f)=b−a2Gn(f∘ϕ)G_n^{[a,b]}(f)=\frac{b-a}2G_n(f\circ\phi)であるから、

∫abf(x) dx−Gn[a,b](f)=b−a2(∫−11f(ϕ(t)) dt−Gn(f∘ϕ))\int_a^bf(x)\,dx-G_n^{[a,b]}(f)=\frac{b-a}2\Bigl(\int_{-1}^1f(\phi(t))\,dt-G_n(f\circ\phi)\Bigr)

である。q∈P2n−1q\in\mathcal P_{2n-1}ならばq∘ϕ∈P2n−1q\circ\phi\in\mathcal P_{2n-1}であるから、(1)により右辺は00である。q(x):=πn(ϕ−1(x))2q(x):=\pi_n(\phi^{-1}(x))^2と置くとq∈P2nq\in\mathcal P_{2n}、q∘ϕ=πn2q\circ\phi=\pi_n^2であり、定理 3.2 (1)により右辺の括弧は∫−11πn(t)2 dt>0\int_{-1}^1\pi_n(t)^2\,dt>0である。したがってGn[a,b]G_n^{[a,b]}の代数的精度は2n−12n-1である。ffがC2nC^{2n}級ならばg:=f∘ϕg:=f\circ\phiは[−1,1][-1,1]上でC2nC^{2n}級であり、g(2n)(t)=(b−a2)2nf(2n)(ϕ(t))g^{(2n)}(t)=\bigl(\frac{b-a}2\bigr)^{2n}f^{(2n)}(\phi(t))である。(1)をggに適用して得るτ∈(−1,1)\tau\in(-1,1)についてξ:=ϕ(τ)∈(a,b)\xi:=\phi(\tau)\in(a,b)と置き、上の等式に代入して主張を得る。▨

例 4.4.n=2n=2とする。P2=(3x2−1)/2P_2=(3x^2-1)/2の零点はx1=−1/3x_1=-1/\sqrt3、x2=1/3x_2=1/\sqrt3であり、L1(x)=32(13−x)L_1(x)=\frac{\sqrt3}2\bigl(\frac1{\sqrt3}-x\bigr)、L2(x)=32(x+13)L_2(x)=\frac{\sqrt3}2\bigl(x+\frac1{\sqrt3}\bigr)からw1=w2=1w_1=w_2=1である。したがってG2(f)=f(−1/3)+f(1/3)G_2(f)=f(-1/\sqrt3)+f(1/\sqrt3)である。k=0,1,2,3k=0,1,2,3に対してG2(xk)G_2(x^k)は2,0,2/3,02,0,2/3,0であり、∫−11xk dx\int_{-1}^1x^k\,dxに等しい。G2(x4)=2/9G_2(x^4)=2/9、∫−11x4 dx=2/5\int_{-1}^1x^4\,dx=2/5であり、その差は8/458/45である。π2=x2−1/3\pi_2=x^2-1/3について∫−11π22 dx=2/5−4/9+2/9=8/45\int_{-1}^1\pi_2^2\,dx=2/5-4/9+2/9=8/45であり、f=x4f=x^4ではf(4)=24f^{(4)}=24であるから、系 4.3 (1)の誤差の式の右辺はξ\xiによらず244!⋅845=845\frac{24}{4!}\cdot\frac8{45}=\frac8{45}である。[0,1][0,1]上のf(x)=exf(x)=e^xについて、G2[0,1](f)=12(e1/2−1/(23)+e1/2+1/(23))G_2^{[0,1]}(f)=\frac12\bigl(e^{1/2-1/(2\sqrt3)}+e^{1/2+1/(2\sqrt3)}\bigr)を5050桁の十進演算で計算すると1.717896378007…1.717896378007\ldotsであり、誤差は∫01f−G2[0,1](f)=3.854504515…×10−4\int_0^1f-G_2^{[0,1]}(f)=3.854504515\ldots\times10^{-4}である。系 4.3 (2)により誤差は(12)5eξ24⋅845=eξ4320\bigl(\frac12\bigr)^5\frac{e^\xi}{24}\cdot\frac8{45}=\frac{e^\xi}{4320}（ξ∈(0,1)\xi\in(0,1)）に等しく、1/4320=2.3148…×10−41/4320=2.3148\ldots\times10^{-4}とe/4320=6.2923…×10−4e/4320=6.2923\ldots\times10^{-4}の間にある。同じffについて、22個の値を用いるT[0,1](f)T_{[0,1]}(f)の誤差は−1.408590857…×10−1-1.408590857\ldots\times10^{-1}、33個の値を用いるS[0,1](f)S_{[0,1]}(f)の誤差は−5.793234175…×10−4-5.793234175\ldots\times10^{-4}である。

例 4.5.n=4n=4とする。§E10.15 定理 4.1により、P0=1P_0=1、P1=xP_1=xから(k+1)Pk+1=(2k+1)xPk−kPk−1(k+1)P_{k+1}=(2k+1)xP_k-kP_{k-1}でP2,P3,P4P_2,P_3,P_4が定まり、P4=(35x4−30x2+3)/8P_4=(35x^4-30x^2+3)/8である。漸化式を微分した(k+1)Pk+1′=(2k+1)(Pk+xPk′)−kPk−1′(k+1)P_{k+1}'=(2k+1)(P_k+xP_k')-kP_{k-1}'により、点xxにおけるP4(x)P_4(x)とP4′(x)P_4'(x)は、多項式を展開せずにxxから計算される。P4P_4は偶関数であり、P4(0)=3/8≠0P_4(0)=3/8\ne0であるから、系 4.1の零点はx1=−x4x_1=-x_4、x2=−x3x_2=-x_3、0<x3<x4<10<x_3<x_4<1を満たす。有理数の演算により

P4(0.33)=2961447160000000>0,P4(0.35)=−4793256000<0,P4(0.86)=−5339310000000<0,P4(0.87)=6888327160000000>0P_4(0.33)=\frac{2961447}{160000000}>0,\quad P_4(0.35)=-\frac{4793}{256000}<0,\quad P_4(0.86)=-\frac{53393}{10000000}<0,\quad P_4(0.87)=\frac{6888327}{160000000}>0

であるから、§D1.12 定理 1.1によりx3∈(0.33,0.35)x_3\in(0.33,0.35)、x4∈(0.86,0.87)x_4\in(0.86,0.87)である。x3x_3に対してJ:=[0.32,0.36]J:=[0.32,0.36]、初期値0.350.35をとる。JJ上でP4′(x)=x(140x2−60)/8<0P_4'(x)=x(140x^2-60)/8<0であり、x(60−140x2)x(60-140x^2)はJJ上で増加するから∣P4′∣≥∣P4′(0.32)∣=1.82656\lvert P_4'\rvert\ge\lvert P_4'(0.32)\rvert=1.82656、また∣P4′′(x)∣=(60−420x2)/8≤∣P4′′(0.32)∣=2.124\lvert P_4''(x)\rvert=(60-420x^2)/8\le\lvert P_4''(0.32)\rvert=2.124である。K:=2.124/(2⋅1.82656)<0.59K:=2.124/(2\cdot1.82656)<0.59と置く。x∈Jx\in Jならば∣x−x3∣<0.03\lvert x-x_3\rvert<0.03であり、§E20.4 補題 4.2により∣NP4(x)−x3∣≤K∣x−x3∣2<0.018∣x−x3∣<0.00054\lvert N_{P_4}(x)-x_3\rvert\le K\lvert x-x_3\rvert^2<0.018\lvert x-x_3\rvert<0.00054であるから、NP4(x)∈JN_{P_4}(x)\in Jである。x4x_4に対してJ:=[0.85,0.88]J:=[0.85,0.88]、初期値0.870.87をとると、同様にJJ上で∣P4′∣≥P4′(0.85)=4.3721875\lvert P_4'\rvert\ge P_4'(0.85)=4.3721875、∣P4′′∣≤P4′′(0.88)=33.156\lvert P_4''\rvert\le P_4''(0.88)=33.156であり、K:=33.156/(2⋅4.3721875)<3.8K:=33.156/(2\cdot4.3721875)<3.8、∣x−x4∣<0.02\lvert x-x_4\rvert<0.02から∣NP4(x)−x4∣<0.076∣x−x4∣<0.00152\lvert N_{P_4}(x)-x_4\rvert<0.076\lvert x-x_4\rvert<0.00152であって、NP4(x)∈JN_{P_4}(x)\in Jである。いずれの場合も、帰納法により Newton 法の反復列(yk)(y_k)はすべてJJに属し、∣yk+1−xi∣≤K∣yk−xi∣2\lvert y_{k+1}-x_i\rvert\le K\lvert y_k-x_i\rvert^2を満たす。P4(x)=0P_4(x)=0はx2x^2の二次方程式35x4−30x2+3=035x^4-30x^2+3=0であるから、x3=(15−230)/35=0.339981043584856…x_3=\sqrt{(15-2\sqrt{30})/35}=0.339981043584856\ldots、x4=(15+230)/35=0.861136311594052…x_4=\sqrt{(15+2\sqrt{30})/35}=0.861136311594052\ldotsである。5050桁の十進演算で計算した反復の誤差yk−xiy_k-x_iは次のとおりである。

kk 00 11 22 33 44
x3x_3（y0=0.35y_0=0.35） 1.002×10−21.002\times10^{-2} 3.188×10−53.188\times10^{-5} 3.904×10−103.904\times10^{-10} 5.858×10−205.858\times10^{-20} 1.319×10−391.319\times10^{-39}
x4x_4（y0=0.87y_0=0.87） 8.864×10−38.864\times10^{-3} 2.512×10−42.512\times10^{-4} 2.100×10−72.100\times10^{-7} 1.470×10−131.470\times10^{-13} 7.199×10−267.199\times10^{-26}

x5−i=−xix_{5-i}=-x_iからLi(−x)=L5−i(x)L_i(-x)=L_{5-i}(x)であり、w1=w4w_1=w_4、w2=w3w_2=w_3である。系 4.3 (1)によりG4G_4は11とx2x^2を厳密に積分するから、2w3+2w4=22w_3+2w_4=2、2w3x32+2w4x42=2/32w_3x_3^2+2w_4x_4^2=2/3であり、

w3=x42−1/3x42−x32=18+3036=0.652145154862546…,w4=18−3036=0.347854845137453…w_3=\frac{x_4^2-1/3}{x_4^2-x_3^2}=\frac{18+\sqrt{30}}{36}=0.652145154862546\ldots,\qquad w_4=\frac{18-\sqrt{30}}{36}=0.347854845137453\ldots

である。

5 Romberg 積分

命題 5.1.a<ba<bを実数、m∈N≥1m\in\NNとし、f ⁣:[a,b]→Rf\colon[a,b]\to\RをC2mC^{2m}級の関数とする。BjB_jを Bernoulli 数（§E5.22 定義 1.1）とし、1≤r≤m1\le r\le mに対して

c2r:=B2r(2r)!(f(2r−1)(b)−f(2r−1)(a)),Dm:=1(2m)!(∣B2m∣+∑j=02m(2mj)∣Bj∣)∫ab∣f(2m)(x)∣ dxc_{2r}:=\frac{B_{2r}}{(2r)!}\bigl(f^{(2r-1)}(b)-f^{(2r-1)}(a)\bigr),\qquad D_m:=\frac{1}{(2m)!}\Bigl(\lvert B_{2m}\rvert+\sum_{j=0}^{2m}\binom{2m}j\lvert B_j\rvert\Bigr)\int_a^b\lvert f^{(2m)}(x)\rvert\,dx

と置く。任意のN∈N≥1N\in\NNに対して、h:=(b−a)/Nh:=(b-a)/Nとすると

∣Th(f)−∫abf(x) dx−∑r=1m−1c2rh2r∣≤Dmh2m\Bigl|T_h(f)-\int_a^bf(x)\,dx-\sum_{r=1}^{m-1}c_{2r}h^{2r}\Bigr|\le D_mh^{2m}

が成り立つ。

証明.g(s):=f(a+hs)g(s):=f(a+hs)（s∈[0,N]s\in[0,N]）と置く。ggはC2mC^{2m}級であり、0≤k≤2m0\le k\le2mに対してg(k)(s)=hkf(k)(a+hs)g^{(k)}(s)=h^kf^{(k)}(a+hs)である。§E5.22 定理 3.1を整数0<N0<Nとggに適用し、C2m:=∑j=02m(2mj)∣Bj∣C_{2m}:=\sum_{j=0}^{2m}\binom{2m}j\lvert B_j\rvertと置くと、

∑j=0Ng(j)=∫0Ng(s) ds+g(0)+g(N)2+∑r=1mB2r(2r)!(g(2r−1)(N)−g(2r−1)(0))+R,∣R∣≤C2m(2m)!∫0N∣g(2m)(s)∣ ds\sum_{j=0}^Ng(j)=\int_0^Ng(s)\,ds+\frac{g(0)+g(N)}2+\sum_{r=1}^m\frac{B_{2r}}{(2r)!}\bigl(g^{(2r-1)}(N)-g^{(2r-1)}(0)\bigr)+R,\qquad\lvert R\rvert\le\frac{C_{2m}}{(2m)!}\int_0^N\lvert g^{(2m)}(s)\rvert\,ds

を満たす実数RRが存在する。置換積分x=a+hsx=a+hsにより∫0Ng(s) ds=h−1∫abf(x) dx\int_0^Ng(s)\,ds=h^{-1}\int_a^bf(x)\,dx、∫0N∣g(2m)(s)∣ ds=h2m−1∫ab∣f(2m)(x)∣ dx\int_0^N\lvert g^{(2m)}(s)\rvert\,ds=h^{2m-1}\int_a^b\lvert f^{(2m)}(x)\rvert\,dxであり、g(2r−1)(N)−g(2r−1)(0)=h2r−1(f(2r−1)(b)−f(2r−1)(a))g^{(2r-1)}(N)-g^{(2r-1)}(0)=h^{2r-1}\bigl(f^{(2r-1)}(b)-f^{(2r-1)}(a)\bigr)である。Th(f)=h(∑j=0Ng(j)−(g(0)+g(N))/2)T_h(f)=h\bigl(\sum_{j=0}^Ng(j)-(g(0)+g(N))/2\bigr)であるから、上の等式にhhを掛けて

Th(f)−∫abf(x) dx−∑r=1m−1c2rh2r=c2mh2m+hR,∣hR∣≤C2m(2m)!h2m∫ab∣f(2m)(x)∣ dxT_h(f)-\int_a^bf(x)\,dx-\sum_{r=1}^{m-1}c_{2r}h^{2r}=c_{2m}h^{2m}+hR,\qquad\lvert hR\rvert\le\frac{C_{2m}}{(2m)!}h^{2m}\int_a^b\lvert f^{(2m)}(x)\rvert\,dx

を得る。§D1.19 定理 2.1によりf(2m−1)(b)−f(2m−1)(a)=∫abf(2m)(x) dxf^{(2m-1)}(b)-f^{(2m-1)}(a)=\int_a^bf^{(2m)}(x)\,dxであるから、§D1.17 定理 3.5により∣c2m∣≤∣B2m∣(2m)!∫ab∣f(2m)(x)∣ dx\lvert c_{2m}\rvert\le\frac{\lvert B_{2m}\rvert}{(2m)!}\int_a^b\lvert f^{(2m)}(x)\rvert\,dxである。二つの評価を加えて主張を得る。▨

定義 5.2.a<ba<bを実数、f ⁣:[a,b]→Rf\colon[a,b]\to\R、N0∈N≥1N_0\in\NNとする。k∈N≥0k\in\Nに対してhk:=(b−a)/(2kN0)h_k:=(b-a)/(2^kN_0)と置く。Thk(f)T_{h_k}(f)は[a,b][a,b]を2kN02^kN_0等分した複合台形則であり、hk+1=hk/2h_{k+1}=h_k/2であるから、kkを11増やすと分割数は22倍になる。Rk,0:=Thk(f)R_{k,0}:=T_{h_k}(f)（k∈N≥0k\in\N）とし、1≤j≤k1\le j\le kに対して

Rk,j:=Rk,j−1+Rk,j−1−Rk−1,j−14j−1R_{k,j}:=R_{k,j-1}+\frac{R_{k,j-1}-R_{k-1,j-1}}{4^j-1}

と置く。(Rk,j)k∈N≥0, 0≤j≤k(R_{k,j})_{k\in\N,\,0\le j\le k}をffの Romberg の表 (Romberg table) といい、Rk,jR_{k,j}（k≥jk\ge j）をその第jj列という。Romberg の表によって∫abf(x) dx\int_a^bf(x)\,dxを近似することを Romberg 積分 (Romberg integration) という。

定理 5.3.a<ba<bを実数、m∈N≥1m\in\NN、N0∈N≥1N_0\in\NNとし、f ⁣:[a,b]→Rf\colon[a,b]\to\RをC2mC^{2m}級の関数、hkh_kとRk,jR_{k,j}を定義 5.2のとおりとする。整数0≤j≤m−10\le j\le m-1に対して、kkによらない定数Kj≥0K_j\ge0が存在して、k≥jk\ge jを満たす任意の整数kkについて

∣Rk,j−∫abf(x) dx∣≤Kjhk2j+2\Bigl|R_{k,j}-\int_a^bf(x)\,dx\Bigr|\le K_jh_k^{2j+2}

が成り立つ。

証明.L:=∫abf(x) dxL:=\int_a^bf(x)\,dx、h∗:=(b−a)/N0h_*:=(b-a)/N_0と置き、c2rc_{2r}とDmD_mを命題 5.1のとおりとする。関数A ⁣:(0,h∗]→RA\colon(0,h_*]\to\Rを、h=(b−a)/Nh=(b-a)/N（N∈N≥1N\in\NN、N≥N0N\ge N_0）のときA(h):=Th(f)A(h):=T_h(f)、それ以外のhhについてA(h):=L+∑r=1m−1c2rh2rA(h):=L+\sum_{r=1}^{m-1}c_{2r}h^{2r}で定める。命題 5.1により、任意の0<h≤h∗0<h\le h_*について∣A(h)−L−∑r=1m−1c2rh2r∣≤Dmh2m\bigl|A(h)-L-\sum_{r=1}^{m-1}c_{2r}h^{2r}\bigr|\le D_mh^{2m}である。すなわちAAは極限LL、係数c2,c4,…,c2m−2c_2,c_4,\dots,c_{2m-2}、定数DmD_mの、指数2,4,…,2m2,4,\dots,2mの誤差の漸近展開をもつ。

m=1m=1ならばj=0j=0であり、Rk,0=A(hk)R_{k,0}=A(h_k)から∣Rk,0−L∣≤D1hk2\lvert R_{k,0}-L\rvert\le D_1h_k^2である。m≥2m\ge2とする。§E20.19 定理 4.4を、係数の個数m−1m-1、指数pr:=2rp_r:=2r（1≤r≤m1\le r\le m）、比t:=2t:=2としてAAに適用する。ここで§E20.19 定理 4.4の係数crc_rはc2rc_{2r}である。A0:=AA_0:=A、Aj:=T2,2jAj−1A_j:=T_{2,2j}A_{j-1}（1≤j≤m−11\le j\le m-1）である。

主張 5.3.1.0≤j≤m−10\le j\le m-1とk≥jk\ge jを満たす整数j,kj,kについてRk,j=Aj(hk−j)R_{k,j}=A_j(h_{k-j})である。

証明.hk=(b−a)/(2kN0)h_k=(b-a)/(2^kN_0)かつ2kN0≥N02^kN_0\ge N_0であるから、Rk,0=Thk(f)=A(hk)R_{k,0}=T_{h_k}(f)=A(h_k)である。1≤j≤m−11\le j\le m-1とし、k′≥j−1k'\ge j-1を満たすすべての整数k′k'についてRk′,j−1=Aj−1(hk′−(j−1))R_{k',j-1}=A_{j-1}(h_{k'-(j-1)})であるとする。k≥jk\ge jならば、§E20.19 定義 4.1の第二の表示とhk−j/2=hk−(j−1)h_{k-j}/2=h_{k-(j-1)}により

Aj(hk−j)=Aj−1(hk−(j−1))+Aj−1(hk−(j−1))−Aj−1(h(k−1)−(j−1))4j−1=Rk,j−1+Rk,j−1−Rk−1,j−14j−1=Rk,jA_j(h_{k-j})=A_{j-1}(h_{k-(j-1)})+\frac{A_{j-1}(h_{k-(j-1)})-A_{j-1}(h_{(k-1)-(j-1)})}{4^j-1}=R_{k,j-1}+\frac{R_{k,j-1}-R_{k-1,j-1}}{4^j-1}=R_{k,j}

である。▨

§E20.19 定理 4.4により、AjA_jは極限LL、係数cj+1(j),…,cm−1(j)c^{(j)}_{j+1},\dots,c^{(j)}_{m-1}、定数CjC_jの、指数2j+2,…,2m2j+2,\dots,2mの誤差の漸近展開をもつ。したがって0<h≤h∗0<h\le h_*について

∣Aj(h)−L∣≤∑r=j+1m−1∣cr(j)∣h2r+Cjh2m≤Kj′h2j+2,Kj′:=∑r=j+1m−1∣cr(j)∣h∗2r−2j−2+Cjh∗2m−2j−2\lvert A_j(h)-L\rvert\le\sum_{r=j+1}^{m-1}\lvert c^{(j)}_r\rvert h^{2r}+C_jh^{2m}\le K'_jh^{2j+2},\qquad K'_j:=\sum_{r=j+1}^{m-1}\lvert c^{(j)}_r\rvert h_*^{2r-2j-2}+C_jh_*^{2m-2j-2}

である。hk−j=2jhk≤h∗h_{k-j}=2^jh_k\le h_*であるから、主張 5.3.1により∣Rk,j−L∣=∣Aj(hk−j)−L∣≤Kj′(2jhk)2j+2\lvert R_{k,j}-L\rvert=\lvert A_j(h_{k-j})-L\rvert\le K'_j(2^jh_k)^{2j+2}であり、Kj:=4j(j+1)Kj′K_j:=4^{j(j+1)}K'_jと置いて主張を得る。▨

例 5.4.f(x):=exf(x):=e^x、[a,b]=[0,1][a,b]=[0,1]、N0=1N_0=1とする。hk=2−kh_k=2^{-k}であり、5050桁の十進演算で計算した∫01f−Rk,j\int_0^1f-R_{k,j}は次のとおりである。

kk j=0j=0 j=1j=1 j=2j=2 j=3j=3 j=4j=4
00 −1.4086×10−1-1.4086\times10^{-1}
11 −3.5649×10−2-3.5649\times10^{-2} −5.7932×10−4-5.7932\times10^{-4}
22 −8.9401×10−3-8.9401\times10^{-3} −3.7013×10−5-3.7013\times10^{-5} −8.5947×10−7-8.5947\times10^{-7}
33 −2.2368×10−3-2.2368\times10^{-3} −2.3262×10−6-2.3262\times10^{-6} −1.3759×10−8-1.3759\times10^{-8} −3.3549×10−10-3.3549\times10^{-10}
44 −5.5930×10−4-5.5930\times10^{-4} −1.4559×10−7-1.4559\times10^{-7} −2.1631×10−10-2.1631\times10^{-10} −1.3434×10−12-1.3434\times10^{-12} −3.3087×10−14-3.3087\times10^{-14}

ffは任意のm∈N≥1m\in\NNについてC2mC^{2m}級であるから、定理 5.3により第jj列の誤差はhk2j+2=4−k(j+1)h_k^{2j+2}=4^{-k(j+1)}の定数倍以下である。第jj列でkkを11増やしたときの誤差の比(∫01f−Rk−1,j)/(∫01f−Rk,j)(\int_0^1f-R_{k-1,j})/(\int_0^1f-R_{k,j})は、j=0j=0で3.951,3.988,3.997,3.9993.951,3.988,3.997,3.999、j=1j=1で15.65,15.91,15.9815.65,15.91,15.98、j=2j=2で62.46,63.6162.46,63.61、j=3j=3で249.7249.7であり、4j+14^{j+1}に近い。k≥1k\ge1とし、h:=hk−1h:=h_{k-1}と置く。4Th/2(f)−Th(f)4T_{h/2}(f)-T_h(f)の各関数値の係数を比べると3Sh/2(f)3S_{h/2}(f)に等しいから、Rk,1=(4Thk(f)−Thk−1(f))/3=Shk(f)R_{k,1}=(4T_{h_k}(f)-T_{h_{k-1}}(f))/3=S_{h_k}(f)である。

6 Chebyshev 補間による求積

定義 6.1.n∈N≥0n\in\Nとし、c0,…,cnc_0,\dots,c_nを[−1,1][-1,1]のn+1n+1個の Chebyshev 節点とする。節点c0,…,cnc_0,\dots,c_nの[−1,1][-1,1]上の補間型求積公式をFnF_nと書き、n+1n+1点の Fejér の第一公式 (Fejér's first rule) という。

命題 6.2.n∈N≥0n\in\Nとし、InI_nを[−1,1][-1,1]のn+1n+1個の Chebyshev 節点における補間作用素、∥⋅∥∞\lVert\cdot\rVert_\inftyを[−1,1][-1,1]上の最大値ノルム、FnF_nをn+1n+1点の Fejér の第一公式とする。

  1. k∈N≥0k\in\Nに対して、∫−11Tk(x) dx\int_{-1}^1T_k(x)\,dxはkkが奇数のとき00、kkが偶数のとき2/(1−k2)2/(1-k^2)である。
  2. f∈C([−1,1])f\in C([-1,1])とし、bk:=ak(Inf)b_k:=a_k(I_nf)を多項式InfI_nfのkk次の Chebyshev 係数とする。このとき Fn(f)=∫−11Inf(x) dx=b0+∑1≤l≤n/22b2l1−4l2F_n(f)=\int_{-1}^1I_nf(x)\,dx=b_0+\sum_{1\le l\le n/2}\frac{2b_{2l}}{1-4l^2} である。
  3. f∈C([−1,1])f\in C([-1,1])とし、bkb_kを(2)のとおりとし、θj:=(2j+1)π/(2n+2)\theta_j:=(2j+1)\pi/(2n+2)（0≤j≤n0\le j\le n）と置く。0≤k≤n0\le k\le nに対して bk=2n+1∑j=0nf(cos⁡θj)cos⁡kθjb_k=\frac2{n+1}\sum_{j=0}^nf(\cos\theta_j)\cos k\theta_j である。
  4. FnF_nの代数的精度はnn以上である。
  5. f∈C([−1,1])f\in C([-1,1])ならば∣∫−11f(x) dx−Fn(f)∣≤2∥f−Inf∥∞\bigl|\int_{-1}^1f(x)\,dx-F_n(f)\bigr|\le2\lVert f-I_nf\rVert_\inftyである。
  6. ρ>1\rho>1とし、U⊆CU\subseteq\Cを Bernstein 楕円Eρ\mathcal E_\rhoを含む開集合、ffをUU上の正則関数でf([−1,1])⊆Rf([-1,1])\subseteq\Rを満たすもの、M:=max⁡w∈∂Eρ∣f(w)∣M:=\max_{w\in\partial\mathcal E_\rho}\lvert f(w)\rvertとする。ffの[−1,1][-1,1]への制限について ∣∫−11f(x) dx−Fn(f)∣≤8Mρ−nρ−1\Bigl|\int_{-1}^1f(x)\,dx-F_n(f)\Bigr|\le\frac{8M\rho^{-n}}{\rho-1} である。

証明.(1)を示す。置換積分x=cos⁡θx=\cos\theta（θ∈[0,π]\theta\in[0,\pi]）と§E20.11 補題 3.2 (2)により

∫−11Tk(x) dx=∫0πcos⁡kθsin⁡θ dθ=12∫0π(sin⁡(1+k)θ+sin⁡(1−k)θ) dθ\int_{-1}^1T_k(x)\,dx=\int_0^\pi\cos k\theta\sin\theta\,d\theta=\frac12\int_0^\pi\bigl(\sin(1+k)\theta+\sin(1-k)\theta\bigr)\,d\theta

である。整数llについて∫0πsin⁡lθ dθ\int_0^\pi\sin l\theta\,d\thetaは、l=0l=0のとき00、l≠0l\ne0のとき(1−(−1)l)/l(1-(-1)^l)/lである。kkが奇数ならば1±k1\pm kは偶数であり、積分は00である。kkが偶数ならば1±k1\pm kは奇数であり、積分は11+k+11−k=21−k2\frac1{1+k}+\frac1{1-k}=\frac2{1-k^2}である。

(2)を示す。InfI_nfは節点c0,…,cnc_0,\dots,c_nにおけるffの補間多項式であるから、命題 1.3 (1)によりFn(f)=∫−11Inf(x) dxF_n(f)=\int_{-1}^1I_nf(x)\,dxである。Inf∈PnI_nf\in\mathcal P_nであるから、§E20.15 命題 1.2 (1)によりInf=b0/2+∑k=1nbkTkI_nf=b_0/2+\sum_{k=1}^nb_kT_kである。(1)によりkkについて積分して等式を得る。

(3)を示す。N:=n+1N:=n+1、φ:=π/(2N)\varphi:=\pi/(2N)と置くとθj=(2j+1)φ\theta_j=(2j+1)\varphiである。整数mmについて、m=0m=0ならば∑j=0ncos⁡mθj=N\sum_{j=0}^n\cos m\theta_j=Nである。0<∣m∣<2N0<\lvert m\rvert<2Nとする。加法定理により2sin⁡(mφ)cos⁡((2j+1)mφ)=sin⁡(2(j+1)mφ)−sin⁡(2jmφ)2\sin(m\varphi)\cos\bigl((2j+1)m\varphi\bigr)=\sin\bigl(2(j+1)m\varphi\bigr)-\sin(2jm\varphi)であるから、j=0,…,nj=0,\dots,nについて和をとると

2sin⁡(mφ)∑j=0ncos⁡mθj=sin⁡(2Nmφ)−sin⁡0=sin⁡mπ=02\sin(m\varphi)\sum_{j=0}^n\cos m\theta_j=\sin(2Nm\varphi)-\sin0=\sin m\pi=0

である。0<∣mφ∣<π0<\lvert m\varphi\rvert<\piであるからsin⁡(mφ)≠0\sin(m\varphi)\ne0であり、∑j=0ncos⁡mθj=0\sum_{j=0}^n\cos m\theta_j=0である。整数0≤k,l≤n0\le k,l\le nをとる。実数θ\thetaについてcos⁡kθcos⁡lθ=12(cos⁡(k+l)θ+cos⁡(k−l)θ)\cos k\theta\cos l\theta=\frac12\bigl(\cos(k+l)\theta+\cos(k-l)\theta\bigr)であり、0≤k+l≤2n<2N0\le k+l\le2n<2N、∣k−l∣≤n<2N\lvert k-l\rvert\le n<2Nである。k+l=0k+l=0はk=l=0k=l=0と同値であり、k−l=0k-l=0はk=lk=lと同値であるから、上の和の値により

∑j=0ncos⁡kθjcos⁡lθj={N(k=l=0)N/2(k=l≥1)0(k≠l)\sum_{j=0}^n\cos k\theta_j\cos l\theta_j=\begin{cases}N&(k=l=0)\\N/2&(k=l\ge1)\\0&(k\ne l)\end{cases}

である。§E20.15 命題 1.2 (1)によりInf=b0/2+∑k=1nbkTkI_nf=b_0/2+\sum_{k=1}^nb_kT_kである。Chebyshev 節点はcj=cos⁡θjc_j=\cos\theta_jであり、Inf(cj)=f(cj)I_nf(c_j)=f(c_j)であるから、§E20.11 補題 3.2 (2)により0≤j≤n0\le j\le nについてf(cos⁡θj)=b0/2+∑k=1nbkcos⁡kθjf(\cos\theta_j)=b_0/2+\sum_{k=1}^nb_k\cos k\theta_jである。0≤l≤n0\le l\le nとし、両辺にcos⁡lθj\cos l\theta_jを掛けてjjについて和をとり、直交関係を用いると、l=0l=0ならば∑jf(cos⁡θj)=Nb0/2\sum_jf(\cos\theta_j)=Nb_0/2、l≥1l\ge1ならば∑jf(cos⁡θj)cos⁡lθj=Nbl/2\sum_jf(\cos\theta_j)\cos l\theta_j=Nb_l/2である。いずれの場合もbl=2N∑j=0nf(cos⁡θj)cos⁡lθjb_l=\frac2N\sum_{j=0}^nf(\cos\theta_j)\cos l\theta_jである。

(4)は命題 1.3 (1)による。

(5)を示す。(2)により∫−11f(x) dx−Fn(f)=∫−11(f(x)−Inf(x)) dx\int_{-1}^1f(x)\,dx-F_n(f)=\int_{-1}^1\bigl(f(x)-I_nf(x)\bigr)\,dxであり、§D1.17 定理 3.5と§D1.17 命題 3.4により右辺の絶対値は2∥f−Inf∥∞2\lVert f-I_nf\rVert_\infty以下である。

(6)は、(5)と§E20.15 定理 6.6 (2)から従う。▨

例 6.3.α>0\alpha>0に対してfα(z):=1/(1+α2z2)f_\alpha(z):=1/(1+\alpha^2z^2)と置く。fαf_\alphaはUα:=C∖{i/α,−i/α}U_\alpha:=\C\setminus\{i/\alpha,-i/\alpha\}上で正則であり、fα([−1,1])⊆Rf_\alpha([-1,1])\subseteq\R、∫−11fα(x) dx=(2/α)arctan⁡α\int_{-1}^1f_\alpha(x)\,dx=(2/\alpha)\arctan\alphaである。ρ>1\rho>1について±i/α∈Eρ\pm i/\alpha\in\mathcal E_\rhoであることはβρ=(ρ−ρ−1)/2≥1/α\beta_\rho=(\rho-\rho^{-1})/2\ge1/\alphaと同値であり、βρ\beta_\rhoはρ\rhoについて狭義単調増加であるから、Eρ⊆Uα\mathcal E_\rho\subseteq U_\alphaであることはρ<ρα:=(1+1+α2)/α\rho<\rho_\alpha:=(1+\sqrt{1+\alpha^2})/\alphaと同値である。したがって命題 6.2 (6)は1<ρ<ρα1<\rho<\rho_\alphaを満たすρ\rhoについて適用される。ρ5=(1+26)/5=1.21980…\rho_5=(1+\sqrt{26})/5=1.21980\ldots、ρ1=1+2=2.41421…\rho_1=1+\sqrt2=2.41421\ldotsである。α=5\alpha=5、nnが奇数のとき、§E20.15 例 6.7により∑k>n∣ak(f5)∣=226⋅ρ5−(n+1)1−ρ5−2\sum_{k>n}\lvert a_k(f_5)\rvert=\frac2{\sqrt{26}}\cdot\frac{\rho_5^{-(n+1)}}{1-\rho_5^{-2}}であり、§E20.15 定理 6.6 (1)と命題 6.2 (5)により∣∫−11f5−Fn(f5)∣\bigl|\int_{-1}^1f_5-F_n(f_5)\bigr|はこの値の44倍以下である。Fn(fα)F_n(f_\alpha)を5050桁の十進演算で計算した誤差の絶対値と、この上界は次のとおりである。

nn α=5\alpha=5の誤差 α=5\alpha=5の上界 α=1\alpha=1の誤差
1111 1.049×10−21.049\times10^{-2} 4.409×10−14.409\times10^{-1} 3.453×10−73.453\times10^{-7}
2121 2.046×10−42.046\times10^{-4} 6.046×10−26.046\times10^{-2} 1.561×10−111.561\times10^{-11}
4141 9.160×10−89.160\times10^{-8} 1.137×10−31.137\times10^{-3} 9.494×10−209.494\times10^{-20}

同じnnで、α=1\alpha=1の誤差はα=5\alpha=5の誤差の10−410^{-4}倍未満である。

7 演習

問題 7.1.a<ba<bを実数とする。[a,b][a,b]上の Riemann 可積分な関数uuであって、すべてのx∈[a,b]x\in[a,b]でu(x)≥0u(x)\ge0、あるx0∈[a,b]x_0\in[a,b]でu(x0)>0u(x_0)>0を満たし、∫abu(x) dx=0\int_a^bu(x)\,dx=0であるものを構成せよ。構成したuuが Riemann 可積分であることと積分の値を示すこと。

解答.

x0∈[a,b]x_0\in[a,b]をとり、u(x0):=1u(x_0):=1、x≠x0x\ne x_0のときu(x):=0u(x):=0と置く。uuは有界であり、非負で、u(x0)>0u(x_0)>0である。n∈N≥1n\in\NNに対して、[a,b][a,b]をnn等分する分割PnP_nをとる。各小区間は長さ(b−a)/n>0(b-a)/n>0をもち、x0x_0と異なる点を含むから、uuの各小区間での下限は00であって、L(u,Pn)=0L(u,P_n)=0である。uuの小区間での上限は、その小区間がx0x_0を含むとき11、含まないとき00である。x0x_0を含む小区間は高々二つであるから、U(u,Pn)≤2(b−a)/nU(u,P_n)\le2(b-a)/nである。ε>0\varepsilon>0に対して2(b−a)/n<ε2(b-a)/n<\varepsilonを満たすnnをとるとU(u,Pn)−L(u,Pn)<εU(u,P_n)-L(u,P_n)<\varepsilonであり、§D1.17 定理 2.4によりuuは Riemann 可積分である。§D1.17 命題 2.3により、任意のnnについて0=L(u,Pn)≤∫abu(x) dx≤U(u,Pn)≤2(b−a)/n0=L(u,P_n)\le\int_a^bu(x)\,dx\le U(u,P_n)\le2(b-a)/nであるから、∫abu(x) dx=0\int_a^bu(x)\,dx=0である。uuはx0x_0で連続でなく、補題 1.2の連続性の仮定を満たさない。▨

問題 7.2.a<ba<bを実数、wwを(a,b)(a,b)上の重み関数、μ0:=∫abw(x) dx\mu_0:=\int_a^bw(x)\,dx、n∈N≥1n\in\NNとし、ziz_i、λi\lambda_i、GnwG_n^wを定義 3.1のとおりとする。δ≥0\delta\ge0とする。[a,b][a,b]上の実数値関数f,f~f,\tilde fがすべての1≤i≤n1\le i\le nで∣f~(zi)−f(zi)∣≤δ\lvert\tilde f(z_i)-f(z_i)\rvert\le\deltaを満たすならば、∣Gnw(f~)−Gnw(f)∣≤μ0δ\lvert G_n^w(\tilde f)-G_n^w(f)\rvert\le\mu_0\deltaであることを示せ。また、[a,b][a,b]上の任意の実数値関数ffに対して、すべての1≤i≤n1\le i\le nで∣f~(zi)−f(zi)∣≤δ\lvert\tilde f(z_i)-f(z_i)\rvert\le\deltaを満たし、∣Gnw(f~)−Gnw(f)∣=μ0δ\lvert G_n^w(\tilde f)-G_n^w(f)\rvert=\mu_0\deltaを満たす[a,b][a,b]上の実数値関数f~\tilde fが存在することを示せ。

解答.

GnwG_n^wの定義によりGnw(f~)−Gnw(f)=∑i=1nλi(f~(zi)−f(zi))G_n^w(\tilde f)-G_n^w(f)=\sum_{i=1}^n\lambda_i\bigl(\tilde f(z_i)-f(z_i)\bigr)である。定理 3.2 (2)により、すべてのiiでλi>0\lambda_i>0であり、∑i=1nλi=μ0\sum_{i=1}^n\lambda_i=\mu_0である。したがって

∣Gnw(f~)−Gnw(f)∣≤∑i=1nλi∣f~(zi)−f(zi)∣≤δ∑i=1nλi=μ0δ\bigl|G_n^w(\tilde f)-G_n^w(f)\bigr|\le\sum_{i=1}^n\lambda_i\bigl|\tilde f(z_i)-f(z_i)\bigr|\le\delta\sum_{i=1}^n\lambda_i=\mu_0\delta

である。ffを[a,b][a,b]上の実数値関数とし、f~(x):=f(x)+δ\tilde f(x):=f(x)+\delta（x∈[a,b]x\in[a,b]）と置く。すべてのiiで∣f~(zi)−f(zi)∣=δ\lvert\tilde f(z_i)-f(z_i)\rvert=\deltaであり、Gnw(f~)−Gnw(f)=δ∑i=1nλi=μ0δG_n^w(\tilde f)-G_n^w(f)=\delta\sum_{i=1}^n\lambda_i=\mu_0\deltaである。▨

前提記事