求積アルゴリズム集(+Excelのワンセルコード)

はじめに
本題
区分求積法
ニュートンコーツ $(Newton Cotes)$ の公式
ロンバーグ $(Romberg)$ 積分
$3$ 次スプライン補間による求積
ガウス・ルジャンドル $(Gauss-Legendre)$ 求積ガウス・クロンロッド $(Gauss-Kronrod)$ 求積
クレンショー・カーティス $(Clenshaw-Curtis)$ 求積
二重指数関数型数値積分公式 $(DE Fomula)$
指数減衰する関数用の半無限区間 $DE$ 公式
$IMT$ 型積分公式
$DE$ 公式と同じ漸近性能を持つ $IMT$ 型積分公式
おわりに
参考文献

応用数学解説

文献あり

求積アルゴリズム集(+Excelのワンセルコード)

数値積分,求積,Excel

324

この著者は初心者として投稿しています。間違いや考慮が足りていない点が含まれている可能性が高いです。見つけたらコメント欄で優しく指摘してあげましょう。

はじめに

求積アルゴリズムをまとめた記事です。随時更新していきます系の記事欲しさに作りました。
式として表示する際の積分範囲は $[a, b]$ とします。
あとExcelのワンセルコードもできるだけ載せておきます。
pythonのコードとかは探せば出てくると思うんですけど、Excelの数式コードってあんまり出てこないので。(pythonより非効率なコードを強いられることもあるので当たり前ではある。)

ワンセルコードのテスト積分は(テキトーに)これにしました。端点に特異性のある若干いやーな定積分です。
$\int_{0}^{π} e^{\cos x} + \sqrt{x} d x \approx 7.68968192506089$
その下には実際に動かした結果と相対誤差を載せています。
$相対誤差 = \frac{| 結果 - 真値 |}{| 真値 |}$
※fx,の後に関数を入力します。色々関数を変えて試してみてください。使い方はExcelまたはGoogleスプレッドシートのセルの数式バーにそのままコピペするだけ。

(求積法には分割回数、次数などを表す $n$ というパラメータがありますが、関数の評価回数とは別物です。詳しく知りたい方は $n$ とか関数を色々変えたりしてそれぞれの性能を調べてみてください。一応関数の評価回数を $n$ で表したものを記しておきます。)

本題

区分求積法

$\int_{a}^{b} f (x) d x \approx \frac{b - a}{n} \sum_{k = 1}^{n} f (a + \frac{(k - 1) \cdot (b - a)}{n})$
求積と言われて大半の方が思いつくのがこれです。通常複合形式(積分範囲を $n$ 個に分けて求積アルゴリズムを適用させた形式)で使われるのでその表示しか載せていません。 $n \to \infty$ にした高校での区分求積法には出てきますが、求積アルゴリズムとしてはあまり使われません。

区分求積法のワンセルコード

=LET(下限,0,上限,PI(),n,50,x,SEQUENCE(n,1,下限,(上限-下限)/n),fx,EXP(COS(x))+SQRT(x),h,(上限-下限)/n,SUM(IFERROR(fx,0)*h))

$n = 50$ に設定 (関数の評価回数: $n = 50$ )
結果: $7.70465739186755$
相対誤差: $1.9475 \times 10^{- 3}$

ニュートンコーツ $(Newton Cotes)$ の公式

$\int_{a}^{b} f (x) \approx \sum_{i = 0}^{m} w_{i} f (a + \frac{i (b - a)}{m}) (単一適用)$
$w_{i} = \int_{a}^{b} \prod_{j = 1, j \neq i}^{m} \frac{m (x - a) - j (b - a)}{(i - j) (b - a)} d x = (- 1)^{m - 1} \frac{_{m} C_{i}}{m!} \int_{0}^{m} \prod_{j = 0, j \neq i}^{m} (t - j) d t$
この公式を(等間隔の) $(m + 1)$ 点ニュートンコーツの公式と呼んだりします。(開型の公式は割愛させていただきます。)
利点： $w_{i}$ をあらかじめ用意しておけば実装が楽であること。
欠点： $m = 8, m \geq 10$ のとき負の重みが現れ桁落ちの原因になる。
　　　高次になると誤差項が大きくなる場合がある。(ルンゲ現象)
$⟶$ $m$ をべらぼうに大きくすればいいわけではない。
(しかし、推奨されないと分かっていてもやってみたいのがサガですよね。実際やってみると単一適用では全く値が合いませんし、複合形式にしても若干影響が出ます。係数表ならここにあるので気になる方はぜひやってみてください。)
複合公式にすれば安定感が増します。

数値積分のはじめの一歩のような存在のニュートンコーツですが、この等間隔に点を取ることをやめて、特に誤差がひどいところの積分区間を分割して繰り返し求積を行う適応型にすればかなりの実用性があります。(後述する紹介するガウス求積とは別物です。)

台形公式

$\begin{array}{rcl} \int_{a}^{b} f (x) d x & \approx & \sum_{k = 1}^{n} [\frac{1}{2} f (a + \frac{(k - 1) \cdot (b - a)}{n}) + \frac{1}{2} f (a + \frac{k (b - a)}{n})] \cdot \frac{b - a}{n} \\ = & [\frac{f (a)}{2} + \sum_{k = 1}^{n - 1} f (a + \frac{k (b - a)}{n}) + \frac{f (b)}{2}] (複合形式) \end{array}$
先ほどの公式で言う $m = 1$ の場合であり、 $w_{0} = w_{1} = \frac{1}{2}$ です。複合形式では端っこ以外の終点と始点が重なるのでこのような式になります。

台形公式のワンセルコード

=LET(下限,0,上限,PI(),n,50,x,SEQUENCE(n+1,1,下限,(上限-下限)/n),fx,EXP(COS(x))+SQRT(x),h,(上限-下限)/n,SUM(MAP(IFERROR(fx,0),SEQUENCE(n+1),LAMBDA(p,q,IF(OR(q=1,q=n+1),p/2,p))))*h)

$n = 50$ に設定(関数の評価回数: $n + 1 = 51$ )
結果: $7.68650060310704$
相対誤差: $4.1371 \times 10^{- 4}$

シンプソンの公式

$\begin{array}{rcl} \int_{a}^{b} f (x) d x & \approx & \sum_{k = 1}^{n} [\frac{1}{3} f (a + \frac{(k - 1) \cdot (b - a)}{n}) + \frac{4}{3} f (a + \frac{(k - \frac{1}{2}) \cdot (b - a)}{n}) + \frac{1}{3} f (\frac{k \cdot (b - a)}{n})] \cdot \frac{b - a}{n} \\ = & [\frac{f (a)}{3} + \sum_{k = 1}^{n} \frac{4}{3} f (a + \frac{(k - \frac{1}{2}) \cdot (b - a)}{n}) + \sum_{k = 1}^{n - 1} \frac{2}{3} f (a + \frac{k \cdot (b - a)}{n}) + \frac{f (b)}{3}] (複合形式) \end{array}$
先ほどの $m = 2$ の場合であり、 $w_{0} = \frac{1}{3}, w_{1} = \frac{4}{3}, w_{2} = \frac{1}{3}$ となります。

シンプソンの公式のワンセルコード

=LET(下限,0,上限,PI(),n,50,x,SEQUENCE(2*n+1,1,下限,(上限-下限)/(2*n)),fx,EXP(COS(x))+SQRT(x),h,(上限-下限)/(2*n),SUM(MAP(IFERROR(fx,0),SEQUENCE(2*n+1),LAMBDA(p,q,IFS(OR(q=1,q=2*n+1),p,MOD(q,2)=0,4*p,MOD(q,2)=1,2*p))))*h/3)

$n = 50$ に設定(関数の評価回数: $2 n + 1 = 101$ )
結果: $7.68922986258012$
相対誤差: $5.8788 \times 10^{- 5}$

シンプソンの $\frac{3}{8}$ 公式

$\begin{array}{rcl} \int_{a}^{b} f (x) d x & \approx & \sum_{k = 1}^{n} [\frac{3}{8} f (a + \frac{(k - 1) \cdot (b - a)}{n}) + \frac{9}{8} f (a + \frac{(k - \frac{2}{3}) \cdot (b - a)}{n}) + \frac{9}{8} f (a + \frac{(k - \frac{1}{3}) \cdot (b - a)}{n}) + \frac{3}{8} f (\frac{k (b - a)}{n})] \\ = & [\frac{3}{8} f (a) + \frac{9}{8} \sum_{k = 1}^{n} f (a + \frac{(k - \frac{2}{3}) \cdot (b - a)}{n}) + \frac{9}{8} \sum_{k = 1}^{n} f (a + \frac{(k - \frac{1}{3}) \cdot (b - a)}{n}) + \frac{6}{8} \sum_{k = 1}^{n - 1} f (a + \frac{k (b - a)}{n}) + \frac{3}{8} f (b)] (複合形式) \end{array}$

先ほどの $m = 3$ の場合であり、 $w_{0} = \frac{3}{8}, w_{1} = w_{2} = \frac{9}{8}, w_{3} = \frac{3}{8}$ となります。

シンプソンの $\frac{3}{8}$ 公式のワンセルコード

=LET(下限,0,上限,PI(),n,50,x,SEQUENCE(3*n+1,1,下限,(上限-下限)/(3*n)),fx,EXP(COS(x))+SQRT(x),h,(上限-下限)/(3*n),SUM(MAP(IFERROR(fx,0),SEQUENCE(3*n+1),LAMBDA(p,q,IFS(OR(q=1,q=3*n+1),p,MOD(q,3)=1,2*p,OR(MOD(q,3)=0,MOD(q,3)=2),3*p))))*3*h/8)

$n = 50$ に設定 (関数の評価回数: $3 n + 1 = 151$ )
結果: $7.68938232170212$
相対誤差: $3.8962 \times 10^{- 5}$

ブールの公式

$\begin{array}{rcl} \int_{a}^{b} f (x) d x & \approx & \sum_{k = 1}^{n} [\frac{14}{45} f (a + \frac{(k - 1) \cdot (b - a)}{n}) + \frac{64}{45} f (a + \frac{(k - \frac{3}{4}) \cdot (b - a)}{n}) + \frac{8}{15} f (a + \frac{(k - \frac{1}{2}) \cdot (b - a)}{n}) + \frac{64}{45} f (a + \frac{(k - \frac{1}{4}) \cdot (b - a)}{n}) + \frac{14}{45} f (a + \frac{k (b - a)}{n})] \\ = & [\frac{14}{45} f (a) + \frac{64}{45} \sum_{k = 1}^{n} f (a + \frac{(k - \frac{3}{4}) \cdot (b - a)}{n}) + \frac{8}{15} \sum_{k = 1}^{n} f (a + \frac{(k - \frac{1}{2}) \cdot (b - a)}{n}) + \frac{64}{45} \sum_{k = 1}^{n} f (a + \frac{(k - \frac{1}{4}) \cdot (b - a)}{n}) + \frac{28}{45} \sum_{k = 1}^{n - 1} f (a + \frac{k (b - a)}{n}) + \frac{14}{64} f (b)] (複合形式) \end{array}$
$m = 4$ の場合。 $w_{0} = w_{4} = \frac{14}{45}, w_{1} = w_{3} = \frac{64}{45}, w_{2} = \frac{8}{15}$

ブールの公式のワンセルコード

=LET(下限,0,上限,PI(),n,50,x,SEQUENCE(4*n+1,1,下限,(上限-下限)/(4*n)),fx,EXP(COS(x))+SQRT(x),h,(上限-下限)/(4*n),SUM(MAP(IFERROR(fx,0),SEQUENCE(4*n+1),LAMBDA(p,q,IFS(OR(q=1,q=4*n+1),14*p,MOD(q,4)=1,28*p,OR(MOD(q,4)=0,MOD(q,4)=2),64*p,MOD(q,4)=3,24*p))))*h/45)

$n = 50$ に設定 (関数の評価回数: $4 n + 1 = 201$ )
結果: $7.68954157908591$
相対誤差: $1.8251 \times 10^{- 5}$

ロンバーグ $(Romberg)$ 積分

これは先ほどの台形公式にリチャードソンの補外を適用させたものです。
まず、刻み幅を倍々で小さくしながら台形公式を適用し、刻み幅が粗い順に $R_{0, 0}, R_{1, 0}, \dots, R_{n, 0}$ とします。
( $R_{0, 0} = (b - a) \cdot \frac{f (a) + f (b)}{2}$ としましょう。)
これらに対して、次の漸化式を考えます。
$R_{i, j} = R_{i, j - 1} + \frac{R_{i, j - 1} - R_{i - 1, j - 1}}{4^{j} - 1} (1 \leq j \leq i \leq n)$
これをネヴィルのアルゴリズムといいます。これを用いて $R_{n, n}$ を求めます。これの良いところは、台形公式という低次の求積法を用いながら、補外法を適用して精度をより良くしているので、被積分関数が十分なめらかであれば、少ない評価回数かつ十分な精度で求められるところです。
早速Excelのコードを...といきたいところですが、Excelのワンセルコードにおいて、三項間以上の漸化式を一般の $n$ に対して解く方法が本当に思いつきません。(恐らくないと思われますが、もしあればご教授願います。)　よって、このロンバーグ積分をワンセルで記述することはできないと思っていたのですが、何とかなったのでここに記しておきます。実用性などはほとんどない(普通のプログラミング言語であれば配列の再定義で簡単に実装できる)ので、導出は省き、結果とコードだけ記しておきます。
$R_{n, n} = \frac{1}{Q_{n}} \sum_{k = 0}^{n} (- 1)^{n - k} \cdot 4^{k (k + 1) / 2} \cdot {(\binom{n}{k})}_{4} R_{k, 0}$
ただし、 $Q_{n} = \prod_{k = 1}^{n} (4^{k} - 1)$
${(\binom{n}{k})}_{q} = \frac{(1 - q^{n}) (1 - q^{n - 1}) \dots (1 - q^{n - k + 1})}{(1 - q) (1 - q^{2}) \dots (1 - q^{k})}$
※この ${(\binom{n}{k})}_{q}$ のことをガウス二項係数、またはq-二項係数といいます。

ロンバーグ積分のワンセルコード

=LET(a,0,b,PI(),n,8,m,2^n,z,SEQUENCE(m-1,1,a+(b-a)/m,(b-a)/m),fx,LAMBDA(x,EXP(COS(x))+SQRT(x)),qbcn,LAMBDA(i,r,q,IF(OR(r=0,r=i),1,PRODUCT(MAKEARRAY(r,1,LAMBDA(v,w,(1-q^(i-r+v))/(1-q^v)))))),C,fx(z),mSeq,SEQUENCE(m-1,1),T,MAKEARRAY(n+1,1,LAMBDA(p,q,(b-a)/2^(p-1)*SUM((fx(a)+fx(b))/2,IF(MOD(mSeq,2^(n-p+1))=0,INDEX(C,mSeq),0)))),qn,PRODUCT(1-4^SEQUENCE(n)),Rn,MAKEARRAY(ROUND((n+1)/2,0),1,LAMBDA(e,f,qbcn(n,e-1,4))),fpow,SCAN(1,SEQUENCE(n+1,1,0),LAMBDA(f,g,f*4^g))/qn,S,MAKEARRAY(n+1,1,LAMBDA(p,q,INDEX(T,p)*INDEX(Rn,INT((n+2)/2-ABS((n+2)/2-p)))*INDEX(fpow,p)*(-1)^(p-1))),SUM(S))

$n = 8$ に設定:　(関数の評価回数: $2^{n} + 1 = 257$ )
結果: $7.68958872044918$
相対誤差: $1.2120 \times 10^{- 5}$

$3$ 次スプライン補間による求積

$\begin{array}{rcl} \int_{x_{1}}^{x_{n}} f (x) d x & \approx & \sum_{i = 1}^{n - 1} \int_{x_{i}}^{x_{i + 1}} f_{i} (x) d x \\ = & \sum_{i = 1}^{n - 1} (x_{i + 1} - x_{i}) (y_{i} + \frac{A_{i}}{2} + \frac{B_{i}}{3} + \frac{C_{i}}{4}) \end{array}$
これは上記の求積法とは扱いが異なります。上記は全て等間隔の点のみで関数を評価していましたが、この求積法は等間隔でない、
いわゆるデータで取れた $(x, y)$ に対して用いるものです。 $x$ を昇順で並べたときに、その順で全ての点を通るように補間を行います。一応何次でも考えられはするのですが先ほどのルンゲ現象のため、最適なのは $3$ 次であると相場が決まっているようです。

3次スプライン補間のルールとして次のようなものがあります。
・点と点の間を1つの三次関数で補間する。
・あるデータ点で共有する2つの三次関数は、その点での1次微分係数と2次微分係数が両者で等しい。

とりあえずここまでの条件で立てられる式はこんな感じです。
$1 \leq i \leq n$ において、 $(x_{i}, y_{i})$ と置き、
$x_{i}$ から $x_{i + 1}$ までを補間する三次関数を
$f_{i} (x) = y_{i} + A_{i} (\frac{x - x_{i}}{x_{i + 1} - x_{i}}) + B_{i} {(\frac{x - x_{i}}{x_{i + 1} - x_{i}})}^{2} + C_{i} {(\frac{x - x_{i}}{x_{i + 1} - x_{i}})}^{3}$
とすると、
$f_{i} (x_{i + 1}) = f_{i + 1} (x_{x + 1})$ から、
$y_{i} + A_{i} + B_{i} + C_{i} = y_{i + 1}$

$f_{i}^{'} (x_{i + 1}) = f_{i + 1}^{'} (x_{i + 1})$ から、
$\frac{1}{x_{i + 1} - x_{i}} (A_{i} + 2 B_{i} + 3 C_{i}) = \frac{A_{i + 1}}{x_{i + 2} - x_{i + 1}}$

$f_{i}^{″} (x_{i + 1}) = f_{i + 1}^{″} (x_{i + 1})$ から、
$\frac{1}{(x_{i + 1} - x_{i})^{2}} (2 B_{i} + 6 C_{i}) = \frac{B_{i + 1}}{(x_{i + 2} - x_{i + 1})^{2}}$

しかしこれだけですと係数に対して式が2つ足りないので、その2つを何らかの形で補う必要があります。
今回はその補い方を2つ紹介し、そのExcelのコードを記述します。

端点ピン止めスプライン $(Clamped Spline)$

公式の日本語訳が分からないので適当に訳してます。
端点の一次微分係数を指定してあげます。始点の一次微分係数をdif_a、終点の一次微分係数をdif_bとします。
dif_a $= \frac{A_{1}}{x_{2} - x_{1}}$
dif_b $= \frac{1}{x_{n} - x_{n - 1}} (A_{n - 1} + 2 B_{n - 2} + 3 C_{n - 1})$
で表されるので、あとはこれらの式を用いて係数を導出すれば、積分値が求まります。(今回は特にテスト関数による例示をしていません。)
具体的なやり方は秋元幸生「プログラミング言語Pascalとその応用」（工学図書, 1984年）pp.133-140やこちらやこちらが分かりやすいものと思われます。

3次端点ピン止めスプライン補間による求積のワンセルコード

=LET(x,[ここに $x$ のデータを縦一列で入力],y,[ここに $y$ のデータを縦一列で入力],dif_a,[始点の一次微分係数],dif_b,[終点の一次微分係数],n,COUNT(x),dx,MAKEARRAY(n-1,1,LAMBDA(a,b,INDEX(x,a+1)-INDEX(x,a))),dy,MAKEARRAY(n-1,1,LAMBDA(a,b,INDEX(y,a+1)-INDEX(y,a))),xd,MAKEARRAY(n-2,1,LAMBDA(a,b,INDEX(dx,a)/INDEX(dx,a+1))),num,MAKEARRAY(n-2,3,LAMBDA(a,b,COMPLEX(a,b))),k,SCAN(1,num,LAMBDA(a,b,IFS(IMAGINARY(b)=1,(-1)*INDEX(xd,IMREAL(b))/(2-a),IMAGINARY(b)=2,(-1)*INDEX(xd,IMREAL(b))^2/(INDEX(xd,IMREAL(b))-a),IMAGINARY(b)=3,1/(1-a)))),l,SCAN(INDEX(dy,1)-dif_a*INDEX(dx,1),num,LAMBDA(a,b,IFS(IMAGINARY(b)=1,IF(IMREAL(b)=1,-(INDEX(dy,IMREAL(b))+a),-(INDEX(dy,IMREAL(b))+a)/(2-INDEX(k,IMREAL(b)-1,3))),IMAGINARY(b)=2,(INDEX(dy,IMREAL(b))-a)/(INDEX(dx,IMREAL(b))/INDEX(dx,IMREAL(b)+1)-INDEX(k,IMREAL(b),1)),IMAGINARY(b)=3,(INDEX(dy,IMREAL(b)+1)-a)/(1-INDEX(k,IMREAL(b),2))))),M,SCAN((dif_b*INDEX(dx,n-1)-INDEX(dy,n-1)-INDEX(l,n-2,3))/(2-INDEX(k,n-2,3)),num,LAMBDA(a,b,IFS(IMAGINARY(b)=1,IF(IMREAL(b)=1,a,INDEX(l,n-IMREAL(b),1)-INDEX(k,n-IMREAL(b),1)*a),IMAGINARY(b)=2,INDEX(l,n-1-IMREAL(b),3)-INDEX(k,n-1-IMREAL(b),3)*a,IMAGINARY(b)=3,INDEX(l,n-1-IMREAL(b),2)-INDEX(k,n-1-IMREAL(b),2)*a))),cone,INDEX(l,1,1)-INDEX(k,1,1)*INDEX(M,n-2,3),aone,dif_a*(INDEX(dx,1)),bone,INDEX(dy,1)-cone-aone,ycba,MAKEARRAY(n-1,4,LAMBDA(p,q,IF(q=1,INDEX(dx,p)*INDEX(y,p),1/q*INDEX(dx,p)*INDEX(VSTACK(M,HSTACK(cone,bone,aone)),n-p,5-q)))),S,SUM(ycba),S)

自然スプライン $(Natural Spline)$

自然スプラインは始点と終点の二次微分係数を $0$ にとるようなスプライン補間です。一次微分係数を入力する必要がないので楽です。

3次自然スプライン補間による求積のワンセルコード

=LET(x,[ここに $x$ のデータを縦一列で入力],y,[ここに $y$ のデータを縦一列で入力],n,COUNT(x),dx,MAKEARRAY(n-1,1,LAMBDA(a,b,INDEX(x,a+1)-INDEX(x,a))),dy,MAKEARRAY(n-1,1,LAMBDA(a,b,INDEX(y,a+1)-INDEX(y,a))),xd,MAKEARRAY(n-2,1,LAMBDA(a,b,INDEX(dx,a)/INDEX(dx,a+1))),num,MAKEARRAY(n-2,3,LAMBDA(a,b,COMPLEX(a,b))),k,SCAN(0,num,LAMBDA(a,b,IFS(IMAGINARY(b)=1,(-1)*INDEX(xd,IMREAL(b))/(2-a),IMAGINARY(b)=2,(-1)*INDEX(xd,IMREAL(b))^2/(INDEX(xd,IMREAL(b))-a),IMAGINARY(b)=3,1/(1-a)))),l,SCAN(0,num,LAMBDA(a,b,IFS(IMAGINARY(b)=1,IF(IMREAL(b)=1,-(INDEX(dy,IMREAL(b))+a)/2,-(INDEX(dy,IMREAL(b))+a)/(2-INDEX(k,IMREAL(b)-1,3))),IMAGINARY(b)=2,(INDEX(dy,IMREAL(b))-a)/(INDEX(dx,IMREAL(b))/INDEX(dx,IMREAL(b)+1)-INDEX(k,IMREAL(b),1)),IMAGINARY(b)=3,(INDEX(dy,IMREAL(b)+1)-a)/(1-INDEX(k,IMREAL(b),2))))),M,SCAN(INDEX(l,n-2,3)/(INDEX(k,n-2,3)-3),num,LAMBDA(a,b,IFS(IMAGINARY(b)=1,IF(IMREAL(b)=1,a,INDEX(l,n-IMREAL(b),1)-INDEX(k,n-IMREAL(b),1)*a),IMAGINARY(b)=2,INDEX(l,n-1-IMREAL(b),3)-INDEX(k,n-1-IMREAL(b),3)*a,IMAGINARY(b)=3,INDEX(l,n-1-IMREAL(b),2)-INDEX(k,n-1-IMREAL(b),2)*a))),cone,INDEX(l,1,1)-INDEX(k,1,1)*INDEX(M,n-2,3),aone,INDEX(dy,1)-cone,ycba,MAKEARRAY(n-1,4,LAMBDA(p,q,IF(q=1,INDEX(dx,p)*INDEX(y,p),1/q*INDEX(dx,p)*INDEX(VSTACK(M,HSTACK(cone,0,aone)),n-p,5-q)))),S,SUM(ycba),S)

ガウス・ルジャンドル $(Gauss-Legendre)$ 求積ガウス・クロンロッド $(Gauss-Kronrod)$ 求積

まずはガウス・ルジャンドル求積の方から
$\int_{a}^{b} f (x) d x \approx \frac{b - a}{2} \sum_{i = 1}^{m} w_{i} f (\frac{b - a}{2} x_{i} + \frac{a + b}{2}) (単一適用)$
$P_{m} (x)$ を $m$ 次ルジャンドル多項式としたときに、
$x_{i} = i$ 番目に小さい $P_{m} (x) = 0$ の解
$w_{i} = \frac{2 (1 - x_{i}^{2})}{n P_{m - 1} (x_{i})}$

先ほどのニュートン・コーツの公式の方では、 $m$ 次の式に対して高々 $m + 1$ 次 $[m = 偶数の場合]$ の多項式近似しかできませんでした。今回のガウス・ルジャンドル求積は、 $m$ 個の点で $2 m - 1$ 次の多項式まで厳密に積分ができるという優れモノなのです。
この求積法にはルジャンドル多項式という直交多項式が関わっています。こちらが一番詳しく、また分かりやすい説明となっています。直交多項式の性質の良さがばっちり現れています。

注意点として、ニュートンコーツの公式と同じように、あまりに高次のものを使いすぎるとルンゲ現象が起きるかもしれないので $\dots$
っていう話をしようと思ったらこんな論文を見つけてしまいました。振動する特定の関数では、超高次のガウス・ルジャンドル求積を使った方が遥かに効率的であることを示す論文です。なかなかおもしろいですね。
といってもこれは、一部の関数についてであるので、あくまでその使い分けは経験と勘に任されているとしか言えない状況でもあります。

しかしながらこのガウス・ルジャンドル求積、誤差評価がかなりしづらいのが難点です。余剰項から考える手法を取るにしても、関数の $2 n$ 階微分が必要で、実用的とは言えません。かといって次数の違うガウス・ルジャンドル求積を使うと、関数の評価点が異なるため結果を使いまわせず、せっかく少ない評価点で精度よく積分できるという良さがかなり損なわれてしまうことになります。

関数の評価点を使いまわしつつ、誤差評価ができる求積法として、ガウス・クロンロッド求積が考えられました。この求積の分点及び重みに使われている直交多項式はスティルチェイス多項式です。
正直、この求積法を単体で用いる利点は、ガウス・ルジャンドル求積の前ではあまりないと思います。ガウス・クロンロッド求積の $2 m + 1$ 点で $3 m + 1$ 次の多項式まで厳密に積分できますが、これはガウス・ルジャンドル求積の性能には劣ります。だから、ガウス・ルジャンドル求積の補助として使われることが多いわけです。

言い忘れましたがこの求積法の分点と重みの計算が重いので、あらかじめ計算されたものを使ったりすることが多いです。pythonとか julia とかだと専用のパッケージとかあるかもしれませんね。

では、ワンセルコードを示します。今回はよく使われている $7$ 点ルジャンドル/ $15$ 点クロンロッド求積で。
その他の次数の分点及び重みが知りたい方はこちらへ。

$7$ 点ガウス・ルジャンドル求積のワンセルコード

=LET(a,0,b,PI(),n,20,xl,VSTACK(-0.949107912342759,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342759),x,MAKEARRAY(n,7,LAMBDA(s,t,a+(b-a)/(2*n)*(INDEX(xl,t)+2*s-1))),fx,EXP(COS(x))+SQRT(x),sfx,BYCOL(IFERROR(fx,0),LAMBDA(u,SUM(u)))*(b-a)/(2*n),wl,HSTACK(0.12948496616887,0.279705391489277,0.381830050505119,0.417959183673469,0.381830050505119,0.279705391489277,0.12948496616887),gl,SUM(MAP(sfx,wl,LAMBDA(p,q,p*q))),gl)

$n = 20$ に設定 (関数の評価回数: $7 n = 140$ )
結果: $7.68969726603681$
相対誤差: $1.9950 \times 10^{- 6}$

$15$ 点ガウス・クロンロッド求積 (誤差評価付)のワンセルコード

=LET(下限,0,上限,PI(),n,10,xk,VSTACK(-0.991455371120813,-0.949107912342759,-0.864864423359769,-0.741531185599394,-0.586087235469691,-0.405845151377397,-0.207784955007898,0,0.207784955007898,0.405845151377397,0.586087235469691,0.741531185599394,0.864864423359769,0.949107912342759,0.991455371120813),x,MAKEARRAY(n,15,LAMBDA(s,t,下限+(上限-下限)/(2*n)*(INDEX(xk,t)+2*s-1))),fx,EXP(COS(x))+SQRT(x),sfx,BYCOL(IFERROR(fx,0),LAMBDA(u,SUM(u)))*(上限-下限)/(2*n),wk,HSTACK(0.022935322010529,0.063092092629979,0.10479001032225,0.140653259715526,0.169004726639268,0.190350578064785,0.204432940075299,0.209482141084728,0.204432940075299,0.190350578064785,0.169004726639268,0.140653259715526,0.10479001032225,0.063092092629979,0.022935322010529),wl,HSTACK(0,0.12948496616887,0,0.279705391489277,0,0.381830050505119,0,0.417959183673469,0,0.381830050505119,0,0.279705391489277,0,0.12948496616887,0),gk,SUM(MAP(sfx,wk,LAMBDA(p,q,p*q))),gl,SUM(MAP(sfx,wl,LAMBDA(p,q,p*q))),CONCAT(gk,"±",(200*ABS(gk-gl))^(1.5)))

※推定誤差は、 $(200 | G L 7 - G K 15 |)^{1.5}$ で与えられる。

$n = 10$ に設定 (関数評価回数: $15 n = 150$ )
結果: $7.68968429498143 \pm 0.000743109748821089$
絶対誤差: $0.00000236992 . . .$
相対誤差: $3.0819 \times 10^{- 7}$

この誤差推定は、基本的に真の誤差よりも大きく出るので、信頼性が持てます。(例外はあるかもしれない)

ガウス求積は、他の直交多項式を用いることが出来るので様々なバリエーションがありますが、現在は割愛させていただきます。また追記するかもしれません。

クレンショー・カーティス $(Clenshaw-Curtis)$ 求積

今回は分点が偶数( $2 n$ )の場合に限らせていただきます。
$\int_{a}^{b} f (x) d x \approx \sum_{k = 0}^{2 n} \frac{b - a}{2} w_{k} f (x_{k})$
$x_{k} = \cos (\frac{k π}{2 n}), w_{k} = w_{2 n - k} = \frac{2}{n} [\frac{1}{2} + \frac{(- 1)^{k}}{2 - 8 n^{2}} + \sum_{i = 1}^{n - 1} \frac{1}{1 - 4 i^{2}} \cos (\frac{i k π}{n})]$
特に、 $w_{0} = w_{2 n} = \frac{1}{4 n^{2} - 1}$
クレンショー・カーティス求積は分点と重みの計算が比較的軽く、また実用上においてガウス・ルジャンドル求積と劣らない性能を持つ求積です。また、分点の数を $n$ から $2 n$ に増やしたとき、その $2 n$ の点の中に元の $n$ 点が全て含まれるという性質を持っているため誤差評価の効率が良く、多次元の積分まで応用が利きます。(今回のワンセルコードには、誤差評価を載せていません。)

クレンショー・カーティス求積のワンセルコード

=LET(a,0,b,PI(),n,75,x,COS(SEQUENCE(2*n+1,1,0)/(2*n)*PI())*(b-a)/2+(a+b)/2,fx,EXP(COS(x))+SQRT(x),w,MAKEARRAY(n+1,1,LAMBDA(s,t,IF(s=1,1/(4*n^2-1),2/n*SUM(MAKEARRAY(1,n+1,LAMBDA(p,q,IFS(q=1,1/2,q=n+1,(-1)^(s-1)/(1-4*n^2)/2,AND(q<>1,q<>n+1),COS((q-1)*(s-1)*PI()/n)/(1-4*(q-1)^2)))))))),u,SUM(MAKEARRAY(2*n+1,1,LAMBDA(g,h,INDEX(IFERROR(fx,0),g)*INDEX(w,n+1-ABS(n+1-g)))))*(b-a)/2,u)

$n = 75$ に設定 (関数の評価回数: $2 n + 1 = 151$ )
結果: $7.68968174577741$
相対誤差: $2.3315 \times 10^{- 8}$

二重指数関数型数値積分公式 $(DE Fomula)$

$\int_{a}^{b} f (x) d x \approx \sum_{i = 1}^{n} f (\frac{b - a}{2} \tanh (\frac{π}{2} \sinh (t_{i})) + \frac{a + b}{2}) \cdot \frac{\frac{π}{2} \cosh (t_{i})}{\cosh^{2} (\frac{π}{2} \sinh (t_{i}))} \cdot \frac{b - a}{2} h$
ただし、 $t_{i} = - t_{a} + (i - 1) h, h = \frac{2 t_{a}}{n - 1}$
$t_{a}$ は打ち切り範囲とし、倍精度ならば $t_{a}$ は $3.5$ から $4$ くらいが適切だと思われる。

最推しの求積です。分点、重みの計算が軽いのはさることながら、今回の積分のような、端点に特異性のある関数の積分にめっぽう強いという特性を持ちます。あとどんな関数でも比較的精度良く計算できるため汎用性があります。
少し注意点があるとすれば万能ではないことと、分点の数を多くしすぎると累積の誤差が大きくなって若干有効数字が落ちます。
そして端点が(正負関わらず)無限大に発散する場合、特に桁落ちが顕著なので、その場合は変数変換等を行って適切に対処しましょう。
(何の関数かは忘れてしまいましたが、求積の結果が真の結果と大きくずれて収束してしまって、幻影解が現れたことがあります。)

これは有限区間のものですが、半無限区間や無限区間用の公式も存在します。
半無限区間:
下に有界:
$\int_{a}^{\infty} f (x) d x \approx h \sum_{i = 1}^{n} f (e^{\frac{π}{2} \sinh (t_{i})} + a) \cdot e^{\frac{π}{2} \sinh (t_{i})} \frac{π}{2} \cosh (t_{i})$
上に有界:
$\int_{- \infty}^{b} f (x) d x \approx h \sum_{i = 1}^{n} f (- e^{\frac{π}{2} \sinh (t_{i})} + b) \cdot e^{\frac{π}{2} \sinh (t_{i})} \frac{π}{2} \cosh (t_{i})$
無限区間:
$\int_{- \infty}^{\infty} f (x) d x \approx h \sum_{i = 1}^{n} f (\sinh (\frac{π}{2} \sinh (t_{i}))) \cdot \cosh (\frac{π}{2} \sinh (t_{i})) \frac{π}{2} \cosh (t_{i})$
ただし、 $t_{i} = - t_{a} + (i - 1) h$ 、 $h = \frac{2 t_{a}}{n - 1}$

これらもなかなか便利です。
ではワンセルコードを示していきます。

二重指数関数型数値積分公式(有限区間)のワンセルコード

=LET(a,0,b,PI(),n,150,ta,3.5,h,2*ta/(n-1),ti,-ta+(SEQUENCE(n,1,0)*h),x,(b-a)/2*TANH(PI()/2*SINH(ti))+(a+b)/2,fx,EXP(COS(x))+SQRT(x),w,PI()/2*COSH(ti)/(COSH(PI()/2*SINH(ti)))^2,S,SUM(MAP(fx,w,LAMBDA(p,q,p*q))),S*(b-a)/2*h)
$t_{a} = 3.5$ に設定
$n = 150$ に設定 (関数の評価回数: $n = 150$ )
結果: $7.68968192506089$
相対誤差: $< ε$ (マシンイプシロン)

これが二重指数関数型数値積分公式の威力です。
下に有界な半無限区間のテスト積分は
$\int_{0}^{\infty} \frac{(\log x)^{2}}{1 + x^{4}} d x = \frac{3 π^{3}}{32 \sqrt{2}} \approx 2.05544517187372$
とします。

二重指数関数型数値積分公式(下に有界な半無限区間)のワンセルコード

=LET(a,0,n,150,ta,4,h,2*ta/(n-1),ti,-ta+SEQUENCE(n,1,0)*h,y,EXP(PI()/2*SINH(ti)),x,y+a,fx,LN(x)^2/(1+x^4),w,MAP(y,PI()/2*COSH(ti),LAMBDA(p,q,p*q)),S,SUM(MAP(IFERROR(fx,0),w,LAMBDA(p,q,p*q))),S*h)

$t_{a} = 4$ に設定
$n = 150$ に設定 (関数の評価回数: $n = 150$ )
結果: $2.05544517187372$
相対誤差: $< ε$ (マシンイプシロン)

上に有界な半無限区間のテスト積分はこれにします。
$\int_{- \infty}^{0} e^{x} \sin (x) d x = - \frac{1}{2}$

二重指数関数型数値積分公式(上に有界な半無限区間)のワンセルコード

=LET(b,0,n,150,ta,4,h,2*ta/(n-1),ti,-ta+SEQUENCE(n,1,0)*h,y,EXP(PI()/2*SINH(ti)),x,-y+b,fx,SIN(x)*EXP(x),w,MAP(y,PI()/2*COSH(ti),LAMBDA(p,q,p*q)),S,SUM(MAP(IFERROR(fx,0),w,LAMBDA(p,q,p*q))),S*h)

$t_{a} = 4$ に設定
$n = 150$ に設定 (関数の評価回数: $n = 150$ )
結果: $- 0.499999999998908$
相対誤差: $2.1840 \times 10^{- 12}$

無限区間のテスト積分はこれにします。
$\int_{- \infty}^{\infty} \frac{d x}{1 + x^{2}} = π \approx 3.14159265358979$

二重指数関数型数値積分公式(無限区間)のワンセルコード

=LET(n,150,ta,4,h,2*ta/(n-1),ti,-ta+SEQUENCE(n,1,0)*h,x,SINH(PI()/2*SINH(ti)),fx,1/(1+x^2),w,COSH(PI()/2*SINH(ti))*PI()/2*COSH(ti),S,SUM(MAP(IFERROR(fx,0),w,LAMBDA(p,q,p*q))),S*h)

$t_{a} = 4$ に設定
$n = 150$ に設定 (関数の評価回数: $n = 150$ )
結果: $3.14159265358979$
相対誤差: $< ε$ (マシンイプシロン)

ここで、上に有界な半無限区間の積分だけ精度が落ちているのが気になりますが、これはDE公式がもたらす離散化誤差によります。被積分関数がそもそも指数関数オーダーで減衰する場合、二重指数的な変数変換を施すと関数の値が0に近いところを多く評価してしまって、肝心の値がある範囲の評価がおろそかになるため若干精度が落ちます。こういった場合、被積分関数が指数減衰するとき専用のDE公式や、一重指数関数型の数値積分公式(SE公式)を用います。
これと同様に、被積分関数が二重指数減衰する場合は台形則で数値積分を行ったり、これまた専用のDE公式があったりします。

指数減衰する関数用の半無限区間 $DE$ 公式

下に有界:
$\int_{a}^{\infty} f (x) d x \approx h \sum_{i = 1}^{n} f (\exp [- t_{i} + \exp (- t_{i})] + a) \cdot (e^{t_{i}} + 1) \exp [- \exp (- t_{i})]$
上に有界:
$\int_{- \infty}^{b} f (x) d x \approx h \sum_{i = 1}^{n} f (- \exp [- t_{i} + \exp (- t_{i})] + b) \cdot (e^{t_{i}} + 1) \exp [- \exp (- t_{i})]$
ただし、 $t_{i} = - t_{a} + (i - 1) h$ , $h = \frac{2 t_{a}}{n - 1}$

指数減衰・下に有界な半無限区間DE公式のワンセルコード

=LET(a,0,n,150,ta,4,h,2*ta/(n-1),ti,-ta+SEQUENCE(n,1,0)*h,x,a+EXP(ti-EXP(-ti)),fx,SIN(x)*EXP(-x),w,(EXP(ti)+1)*EXP(-EXP(-ti)),S,SUM(MAP(IFERROR(fx,0),w,LAMBDA(p,q,p*q))),S*h)

テスト積分: $\int_{0}^{\infty} e^{- x} \sin x d x = \frac{1}{2}$

$t_{a} = 4$ に設定
$n = 150$ に設定 (関数の評価回数: $n = 150$ )
結果: $0.500000000000000$
相対誤差: $< ε$

指数減衰・上に有界な半無限区間DE公式のワンセルコード

=LET(b,0,n,150,ta,4,h,2*ta/(n-1),ti,-ta+SEQUENCE(n,1,0)*h,x,b-EXP(ti-EXP(-ti)),fx,SIN(x)*EXP(x),w,(EXP(ti)+1)*EXP(-EXP(-ti)),S,SUM(MAP(IFERROR(fx,0),w,LAMBDA(p,q,p*q))),S*h)

テスト積分: $\int_{- \infty}^{0} e^{x} \sin x d x = - \frac{1}{2}$

$t_{a} = 4$ に設定
$n = 150$ に設定 (関数の評価回数: $n = 150$ )
結果: $- 0.500000000000000$
相対誤差: $< ε$

$IMT$ 型積分公式

$\int_{a}^{b} f (x) d x \approx (b - a) \cdot h \sum_{k = 1}^{n} f {(b - a) \cdot ϕ (k h) + a} ϕ^{'} (k h)$
ただし、 $h = \frac{1}{n + 1}$ , $ϕ (t) = \frac{1}{Q} \int_{0}^{t} \exp (- \frac{1}{s} - \frac{1}{1 - s}) d s$ , $Q = \int_{0}^{1} \exp (- \frac{1}{s} - \frac{1}{1 - s}) d s$
これがIMT型積分公式です。IMTという名前はこの公式の開発に関わった方々の頭文字(Iri, Moriguchi, Takasawa)を取ったものです。詳細はこちらに載っています。 $ϕ (t)$ は初等関数では表せないので、漸近展開などを考えなければならないのですが、その展開や計算方法までしっかり載っています。しかしながら、これをExcelのワンセルコードで表すのが今のところ出来ないので...(出来たら追記するかもしれませんが)

$DE$ 公式と同じ漸近性能を持つ $IMT$ 型積分公式

$\int_{a}^{b} f (x) d x \approx \frac{b - a}{2} h \sum_{i = 1}^{n} f (ϕ_{m, k} (t_{i}) \frac{b - a}{2} + \frac{a + b}{2}) \cdot ϕ_{m, k}^{'} (t)$
$h = \frac{2}{n + 1}$ , $t_{i} = - 1 + i h$ , $ϕ_{m, k} (t) = \erf (\frac{k}{(1 - t)^{m}} - \frac{k}{(1 + t)^{m}})$ ,
$m = \frac{1}{2} \log (n + 1), k は定数(パラメータ), \erf (t) = \frac{2}{\sqrt{π}} \int_{0}^{t} e^{- s^{2}} d s$
$ϕ^{'} (t) = \frac{2}{\sqrt{π}} (\frac{m k}{(1 - t)^{m + 1}} + \frac{m k}{(1 + t)^{m + 1}}) \exp (- {(\frac{k}{(1 - t)^{m}} - \frac{k}{(1 + t)^{m}})}^{2})$
こちらを実装していきます。Excelには誤差関数 $\erf (t)$ を計算する $ERF ()$ という関数があるのでそれがそのまま使えるのがうれしいところです。しかも、先述のDE公式と同じ漸近性能
$絶対誤差 = \exp (- c \frac{n}{\log n})$
をもちます。
さらに、Excelでは単精度から倍精度で計算することがほとんどなのですが、その範囲ではDE公式より高性能であることが、いくつかの数値計算例によって確認されています。

$DE$ と同じ漸近性能を持つ $IMT$ 型積分公式のワンセルコード

=LET(a,0,b,PI(),n,150,k,2.2,h,2/(n+1),t,-1+SEQUENCE(n)*h,m,LN(n+1)/2,yp,k/(1-t)^m,yq,k/(1+t)^m,y,yp-yq,x,ERF(y)*(b-a)/2+(a+b)/2,fx,EXP(COS(x))+SQRT(x),w,(yp/(1-t)+yq/(1+t))*EXP(-1*y^2),S,SUM(MAP(IFERROR(fx,0),w,LAMBDA(p,q,p*q))),S/SQRT(PI())*h*(b-a)*m)

$k = 2.2$ に設定
$n = 150$ に設定(関数の評価回数: $n = 150$ )
結果: $7.68968192506089$
相対誤差: $< ε$ (マシンイプシロン)

おわりに

求積アルゴリズムは他にも、モンテカルロ積分などさまざまあります。先ほどのガウス求積も含めて追記していく予定です。
本文、ワンセルコードのミスや、もっと計算量を減らせるなどの提案があればお願いします。
また、これに載っていない求積法などありましたら教えて下さると助かります。
皆さんもExcelやGoogleスプレッドシートを開いて手軽に求積してみてはいかがでしょうか。