従って、台形則の誤差(Euler-Maclaurinの公式)は0に近くなり、効率よく積分を計算することができる。分点数をNとして、計算式は次のようになる。[3])DE公式では、対象とする積分を減衰の速い無限積分に変換してから台形則を適用することにより効率よく積分を計算する。N++N–+1)とすると、計算式は次のようになる。φ(ih)無限積分の計算は、下側を-N–、上側をN+で打ち切ることとした。これらの値は許容誤差限界をもとに適切な値を設定する。1]の積分に変換してからこの公式を適用する。半無限区間の積分において被積分関数がゆるやかに減衰するような場合についても、有限区間の場合と同様の考え方でDE公式を導くことができる。文献[3]に示されているプログラムDEHINTでは、被積分関数により3種類の公式を使い分けるようになっている。被積分関数がゆるやかに減衰するような全無限区間の積分についても同様にDE公式を導くことができる。(π/2)cosh(t)cosh((π/2)sinh(t))数値実験まずは、前に実験した例を用いて、変数変換型公式の性能を見てみる。次の積分を計算してみる。47364DE公式ではプログラムの作りによって関数評価回数が変わってくるので、DE公式だけ横軸をXLPackのDefinルーチンを使用したときの実際の関数評価回数とした。IMT公式、DE公式ともに特別に高性能ではなく、Gauss公式には及ばないようである。変数変換型公式としてはIMT公式が最初に提案され、DE公式は後で開発された。DE公式の性能が良いので、現在では変数変換型公式としては主にDE公式が用いられている。次に、同様に前に実験した例を用いて、端点に特異点がある場合について見てみる。π/4x=1が特異点である。前にも見たように補間型公式ではどれもまともな結果が得られないが、IMT公式とDE公式は特異点の影響を受けていないのがわかる。このように、IMT公式とDE公式の特徴は、端点に特異点がある(端点で微係数が定義できない)場合にも適用できることである。積分区間内に特異点がある場合には良い結果が得られない。基本公式」で見たように、関数評価回数が同程度であれば種々の公式の中でGauss公式が最も精度が良かった。計算するのであればこれで良いが、実務的には積分値と同時に推定誤差も知りたい。後で示す自動積分のベースとして使う場合には誤差推定ができなければならない。Gauss公式の弱点は誤差の推定がしにくいことである。計算結果の誤差を推定するためには、例えば、N分点で計算したものと2N分点で計算したものを比較するなどすればよい。次の例(滑らかな関数)で、台形則の計算例を見てみる。0.000147967誤差推定値はRomberg積分の時と同じように次の式から計算した。実際の誤差と推定値はよく一致している。Iprev)/3ここでは、分点数を倍々にするような順序(分割幅は等間隔で半分)で計算を行った。Iは積分の計算値、Iprevは分点数が半分の時の(つまり1つ前の)計算値である。この場合、分点数を倍にした時にもう一回最初から全部の分点で計算する必要はなく、半分はすでに計算済の値がそのまま使える。ところが、Gauss公式では分点が等間隔でないためそのようなことができないのが弱点となっている。誤差推定機能付きのGauss公式これを解決するために、N点Gauss公式の計算値を生かしたまま2N+1点を使った積分公式がGauss-Kronrod公式である。j=1~N+1)と表される。同じN個の点における係数aiはGauss公式の係数wiとは異なる値を持つ。N点Gauss公式GNは2N-1次多項式まで正しい解を与える(2N-1次である)が、Kronrod公式K2N+1は3N+1次となる。同じ関数評価回数で計算できるG2N+1は4N+1次となるので、Kronrod公式は精度的にはGauss公式に劣るが、誤差の推定ができる利点の方が大きい場合に有効である。Gauss-Kronrod公式では、2N+1回の関数評価で2N+1次のGNと3N+1次のK2N+1の両方を求めることができるので、この2つを比べることにより誤差を推定することができる。推定誤差は経験的に次式で表される。K2N+1|)1.5この式による推定誤差は安全側に求められるとされている。Gauss-Kronrod公式の重みおよび分点の具体的な値を下の付表に示す(Gauss公式への追加分のみ)。数値実験Gauss-Kronrod公式を用いて次の積分を求めてみる。47364これは前にも例題に使用した性質のよい関数の例であり、計算結果は次のようになった。緑色のデータは前項の式で求めた推定誤差である。各データ点に表示した数字はNの値である。してみるとKronrod公式でも思ったよりも精度が落ちていないように見える。推定誤差はかなり安全側に算出されているのがわかる。この例題では収束が速すぎてうまく比較できないので、もう少し難しい次の積分を計算してみる。ln(5)これも前に例題に使用した関数で、積分区間外であるが近くに特異点を持つ例である。計算回数がかなり増えたが傾向は同じであった。Kronrod公式を使うと、それほど大きなペナルティを払わずに誤差を推定できることが確かめられた。自動積分ルーチンは、積分区間と要求精度を与えると、誤差推定値が要求精度を満たすまで刻み幅を順次小さくして積分を繰り返すものである。適応自動積分ルーチンは、これを一歩進めて、積分区間を部分区間に分割し、誤差推定値が要求された精度を満たすように部分区間の幅を自動的に調節する。被積分関数がなめらかにゆっくり変化するような区間では大きな幅を用い、誤差推定値が大きい部分では細かい区間に分割して刻み幅を小さくして計算することにより、最小の計算量で要求精度の結果を得ようとするものである。(文献[5])のように誤差推定機能を持つGauss-Kronrod公式のような不等間隔分点の公式を使ったものもある。QAG適応自動積分ルーチンの例としてQAGについて簡単に説明する。QAGは、要求精度を満たすまで自動的に積分区間を部分区間に分割しながらGauss-Kronrod公式の固定分点のルーチンを呼び出すことにより自動積分を行う。XLPackには倍精度版(Dqag)が収録されている。このプログラムの動作の概要は次のとおりである。全区間をGauss-Kronrodルーチンにより積分を行い、積分値Sと推定誤差Eを求める。なったら、その時の積分値の総和Sを最終結果として返すQAGではベースのGauss-KronrodルーチンとしてQK15(N=7)~QK61(N=30)のうちユーザーが指定したものを使用し、これを自動的に変更することはしない。区間を分割する際に計算済の関数値を再利用することもしない。シンプルなアルゴリズムではあるがよく機能する。6ここでは、XLPackに収録されているDqagを使用し、目的誤差を10-5以下としてベースとしてQK15を指定して計算した。トータルの関数評価回数は135回であった。急峻な変化のところで分割されているのがわかる。試しに同じ条件でQK61を指定してみると次のようになった。0.5~1)、トータルの関数評価回数は183回になった。変化が激しい場合にはQK15を選択し分割数を多くし、そうでない場合にはQK61のように高精度のものを選び分割数を少なくするとよさそうである。下表に示す(xiはGauss公式のuiと同じ)。