令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催

二重指数関数型数値積分公式の理論と発展
大浦拓哉

1

はじめに

解析的に計算できない積分の値を求めるには数値計算に頼らないといけない場合があります。このとき必要になるのが数値積分公式であり、多くの数学ソフトウェアではあらかじめパッケージ化されて提供されています。本講座では、そのパッケージがどうなっているのかということを主眼に置きます。
まず前半は、古典的な補間型数値積分公式について解説します。特に最も基本的な台形則がある特殊な場合において最適公式になるという事実はあまり知られていない
ので重点的にお話したいと思います。後半では、コンピュータを用いた高速かつ高精度な計算が可能な二重指数関数型数値積分公式（DE 公式）について、その公式の普及と発展・応用についてご紹介したいと思います。

2

補間型数値積分公式

補間型数値積分公式は被積分関数を多項式で補間することで得られる数値積分公式
であり、コンピュータが普及する前から使われてきました。本節ではまず標本点が等間隔であるニュートン・コーツ公式の中で最も基本的な台形則を説明し、そして補間型数値積分公式の中で最も高精度とされるガウス型公式を説明します。

2.1

台形則

区間 [a, b] で連続な関数 f (x) の積分に対する複合台形則（あるいは単に台形則）は区間を N 等分して関数を N 個の直線で補間して N 個の台形の面積で近似すること
で次のように得られます。
(

N
−1
∑
1
1
TN = hf (a) +f (a + nh) + f (b)
2
2
n=1

)

ここに h = (b − a)/N は刻み幅です。この台形則の近似誤差は古くからよく研究されていて f (x) が [a, b] で C 2m+2 級のときオイラー・マクローリン展開（オイラーの和
公式）

TN −

∫ b

f (x) dx =
a

m
∑
B2j
j=1 (2j)!

()

h2j f (2j−1) (b) − f (2j−1) (a) + O(h2m+2 )が成り立ちます。ここに Bn はベルヌーイ数です。ただしこの式は漸近展開であり C ∞級のとき無限和にすると収束するとは限りません。この展開は有限和で打ち切ったときに N → ∞, h → 0 とすると h2m+2 のオーダーで打ち切り項が小さくなることを意
味します。

1

令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催

通常の台形則の近似誤差はオイラー・マクローリン展開を第 1 項で打ち切った項
h (f 0 (b) − f 0 (a))/12 でおおよその評価ができて N → ∞, h → 0 で O(h2 ) = O(N −2 )の近似誤差になります。
次に少しの計算手順を追加することで精度を上げることを考えてみましょう。
2

(

)

N
−1
N
−1
∑
∑
h0
4T2N − TN00=
f (a) + 4
f (a + 2nh0 ) + f (b)SN =f (a + 2nh + h ) + 2
3
3
n=0
n=1

の近似誤差をオイラー・マクローリン展開から求めると 2 次の項が消えて O(h04 ) となることがわかります。ここに h0 = (b − a)/(2N ) です。この近似はシンプソン則と呼ばれるもので、関数を 2 次関数で補間した積分近似と一致します。同様にして(16S2N − SN )/15 で 4 次の項を消し去り…という具合に繰り返して得られる算法がロ
ンバーグ積分法です。また等間隔の区間分割で高次の多項式で補間して得られる算法
が高次のニュートン・コーツ公式です。これらの高次のニュートン・コーツ公式やロンバーグ積分法は計算量や精度に関してのちに挙げる方法よりも劣っているので実用の数値計算で用いられることはあまりありません。

2.2

ガウス型公式

ガウス型公式は∫ baf (x)w(x) dx ≈N∑

wn f (an )n=1という形の近似で a, b, w(x), N が与えられたとき∫ bk

x w(x) dx =aN∑(k = 0, 1, · · · , 2N − 1)

wn aknn=1を満たすように標本点 a1 , a2 , · · · , aN , 重み w1 , w2 , · · · , wN を定めることで得られる公式です。すなわちガウス型公式は関数 f (x) が任意の 2N − 1 次多項式のときに正しい値を与える公式です。高次のニュートン・コーツ公式との違いは標本点が等間隔ではない点で、標本点を不等間隔にしてその分精度を稼ごうという発想で得られるのがガウス型公式です。ガウス型公式の標本点と重みは直交多項式と関係していて、通常の [−1, 1] 区間の積分 (w(x) = 1, a = −1, b = 1) のときは標本点 a1 , a2 , · · · , aN はル
ジャンドル多項式 PN (x) のゼロ点で与えられ、重みは wn = 2(1 − a2n )/(N PN −1 (an ))2 ,(n = 1, 2, · · · , N ) で与えられます。以下に具体的なガウス型公式の積分と直交多項式の関係を示します。
∫ 1

ルジャンドル・ガウス則 :チェビシェフ・ガウス則 :−1∫ 1

ルジャンドル多項式

f (x) dx,1f (x) √
dx,
2
−1
1
−
x
∫
∞

ラゲール・ガウス則 :チェビシェフ多項式

f (x)e−x dx,

ラゲール多項式

f (x)e−x dx,

エルミート多項式0

エルミート・ガウス則 :

∫ ∞
−∞

2

2

令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催

一般にこれらのガウス型公式は N を大きくすると N の指数関数のオーダーで近似誤差が小さくなります。しかし、N を大きくすると直交多項式のすべてのゼロ点を計算する手間が関数 f (x) を計算する手間を大きく上回るので高速高精度計算には向きません。これらのガウス型公式はいくつかの広義積分の計算が可能ですが、のちに述べる変数変換型公式に比べて応用範囲が限定的です。

台形則が高精度となる積分例

3

通常、台形則の近似誤差は刻み幅 h に対して O(h2 ) で低い精度しか出ませんが、ある種の積分に対しては非常に高精度になって近似誤差は O(exp(−C/h)) の形になります（C は被積分関数に依存する正定数）。その積分は解析的周期関数の 1 周期積分と解析的関数の全無限区間積分の二つです。この事実が広く知られ応用されるようになったのはコンピュータが普及しだした頃の 1960 年代あたりからです。この節では、
この 2 種類の積分について議論します。

3.1

解析的周期関数の 1 周期積分

例として第 1 種ベッセル関数の積分表示∫ 1

J0 (α) =cos(α sin πx) dx
0の計算を考えてみましょう。この積分は周期 1 の関数の 1 周期積分であることがわかります。これに台形則を適用すると

TN = hN−1∑cos(sin πnh),

n=0

h=

1
N

のようになります。具体的に α = 1, N = 8, h = 1/8 として計算すると

J0 (α) ≈ 0.7651976865579665514497
TN ≈ 0.7651976865579665528870
で近似誤差は 1.437 · 10−18 となります。この積分を同じ条件の刻み幅 h0 = 1/8 のシンプソン則で計算したときの近似誤差は 6.281 · 10−8 になり、台形則よりも精度が落ちます。
この場合の台形則が極めて高い精度になることはオイラー・マクローリン展開の近似誤差の項が周期性よりすべて消え去ることからある程度説明がつきます。これは誤差が完全に 0 になるという意味ではなくて、O(h2m+2 ) よりも小さいオーダーになることを意味します。より正確な誤差の見積もりは複素関数論を用いて行います。もし
実軸を含む幅 2d の帯状領域 {z; |=z| ≤ d} で f (z) が正則でその境界の積分を
∫ b

A=
a

(|f (x + id)| + |f (x − id)|) dx

3令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催とするとき 1 周期積分の台形則の近似誤差は

TN −

∫ b
a

f (x) dx ≤A
exp(2πd/h) − 1と評価できます。この右辺は O(exp(−2πd/h)) で通常の台形則の近似誤差である O(h2 )と比べて極端に小さいオーダーになることがわかります。さらに解析的周期関数の 1
周期積分の場合、すべての積分則の中で標本点数が無限大の極限で漸近的に台形則が最良となる（すべての積分側の中での近似誤差のオーダーの下限を達成する）ことが知られています [3, 6]。

3.2

解析的関数の全無限区間積分

まず全無限区間積分の台形則を∞∑Ih = hf (nh)

n=−∞で定義します。実際の数値計算では有限和(N )Ih = hN+∑f (nh),

N = N− + N+ + 1

n=−N−で計算することに注意します。このときの打ち切り N− , N+ は適切に選ぶ必要があるのですが、この節では話を単純にするため無限和での近似 Ih のみを考えます。
例としてガウス積分
∫ ∞
√
2
π=
e−x dx−∞を考えてみましょう。これに台形則を適用するとIh = h∞∑e−(nh)

2n=−∞のようになります。具体的に h = 1/2 として計算すると√
π ≈ 1.772453850905516027298
Ih ≈ 1.772453850905516052670
で近似誤差は 2.537 · 10−17 となります。
この場合もオイラー・マクローリン展開の式で h を固定して a → ∞, b → −∞ とすることで誤差の項がすべて消え去ることからある程度説明がつきます。より詳細な近似誤差は前小節同様に複素関数論を用いて解析します。もし実軸を含む幅 2d の帯状領域 {z; |=z| ≤ d} で f (z) が正則でlim

max |f (x + iy)| = 0x→±∞ −d≤y≤d

かつ
0

A =

∫ ∞
−∞

(|f (x + id)| + |f (x − id)|) dx < +∞4

令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催

ならば全無限区間の台形則の誤差は

Ih −

∫ ∞
−∞

f (x) dx ≤A0
exp(2πd/h) − 1と評価できます。この全無限区間の場合もすべての全無限区間積分則の中で台形則が漸近的に最良である（標本点密度が ∞ の極限で近似誤差のオーダーの下限を達成する）ことが知られています。

4

変数変換型数値積分公式

変数変換型数値積分公式とは積分を変数変換により性質の良い積分に変換してそれに既存の積分則を適用することで得られる公式です。狭義的には変数変換で 1 周期積分あるいは全無限区間積分に変換してから台形則を適用することで得られる公式を指します。具体的に x = ϕ(t) で変数変換する変数変換型公式は∫ b∫

f (x) dx =aΩf (ϕ(t))ϕ0 (t) dt ≈ h∑f (ϕ(nh))ϕ0 (nh)nという近似になり、{ϕ(nh)} が変数変換型公式の標本点、{hϕ0 (nh)} が変数変換型公式の重みになります。この変数変換で得られる公式の従来の積分公式にはない最大の特徴は関数 f (x) が積分端点で非有界となる広義積分の計算ができる点です。端点で定符号で発散する関数の積分は適当な変数変換で有界な関数の積分に変換されます（ただし振動しながら発散するような場合は振動積分として別に扱う必要があります）。またこの公式は構造が単純で ϕ(t) が計算しやすい関数ならば高速で計算できる特徴があります。ただこの公式で問題となるのは ϕ(t) の選び方で、それが変数変換
型公式の性能を決定づけます。

4.1

1 周期積分に変換する公式

数値計算したい積分が簡単に解析的周期関数の 1 周期積分に変換できる場合（完全楕円積分など）はそれに台形則を適用するのが最適な方法になります。一般の積分の場合は 1970 年に提案された伊理・森口・高澤による IMT 公式が有効になります[4, 13, 3]。この IMT 公式は (0, 1) 区間の積分を (0, 1) 区間の積分に変換する変数変換（IMT 変換）
(

1 ∫t
1
1
ϕIMT (t) =exp − −
Q 0
u 1−u

)

∫ 1

du,

(

1
1
Q=
exp − −
u 1−u
0

)

du

を用いて得られる数値積分公式です（区間が (0, 1) 以外の有限積分は区間を (0, 1) にする線形の変数変換を組み合わせます）。関数 f (x) が C ∞ 級であれば変数変換された関数 f (ϕIMT (t))ϕ0IMT (t) は積分端点 t = 0, 1 ですべての微係数が 0 となって、C ∞ 級の周
期 1 の周期関数とみなすことができます。しかし f (x) が解析的で複素関数として考えると f (ϕIMT (t))ϕ0IMT (t) は t = 0, 1 で真性特異点となり、実軸を含む帯状領域で解析的とはならず近似誤差のオーダーは少し大きく O(exp(−C/h1/2 )) = O(exp(−CN 1/2 ))の形になります。

5令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催

このオーダーは IMT 変換を繰り返し適用することで簡単に改善できて、IMT-double,
IMT-triple, …（変数変換は ϕIMT (ϕIMT (t)), ϕIMT (ϕIMT (ϕIMT (t))), …）に対する近似
誤差のオーダーはそれぞれ O(exp(−CN/ log2 N )), O(exp(−CN/(log N log2 log N ))),…のようになります [10]。この IMT 変換はさらなる改良型がいくつか知られています [8, 30]。

4.2

全無限区間積分に変換する公式

積分を全無限区間に変換して台形則で計算する方法は 1960 年代以降に様々な人が提案しています [1, 2, 5]。例えば (−1, 1) 区間を (−∞, ∞) 区間にする変換

ϕtanh (t) = tanh t = −1 +2
1 + e−2tで得られる tanh 則（一重指数関数型公式、略して SE 公式）や

2 ∫ t −u2eϕerf (t) = erf t = √duπ 0で得られる erf 則や()π
sinh t2で得られる tanh-sinh 則（二重指数関数型公式、略して DE 公式）などがあります。近似公式は
ϕDE (t) = tanh

∫ 1
−1

f (x) dx =∫ ∞f (ϕ(t))ϕ0 (t) dt = I−∞∞∑

f (ϕ(nh))ϕ0 (nh) = Ih

≈ h

n=−∞

≈ h

N+
∑

f (ϕ(nh))ϕ0 (nh) = Ih

(N )

n=−N−

のようになります（区間が (−1, 1) 以外の有限積分は区間を (−1, 1) にする線形の変数変換を組み合わせます）。この公式で問題となるのは無限和を有限和で打ち切るときの打(N )ち切り誤差 |Ih −Ih | です（この公式の近似誤差には離散化誤差 |I −Ih | と打ち切り誤差(N )|Ih − Ih | の 2 つが含まれます）。打ち切り誤差を小さくするには変数変換を急峻な関数に選べばよいのですが、あまり急峻すぎると解析性が損なわれて離散化誤差 |I − Ih |が大きくなります。この 2 つの誤差がともに最小となるのは変数変換後の関数の減少が二重指数関数的であるとき（t → ±∞ の漸近形が f (ϕ(t))ϕ0 (t) ∼ exp(−c exp(|t|))のとき）になされ、h に対して N− , N+ を最適に選んだときの全近似誤差のオーダー
は O(exp(−CN/ log N )), N = N− + N+ + 1 となります。

5

DE 公式

DE 公式は 1974 年にある種の最適な変数変換型公式として高橋・森により提案されたものです [7, 12, 14]。それ以前にも DE 公式は SE 変換の組み合わせとして計算例
6令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催が示されていますが最適性は示されていませんでした [2]。その基本的な指針は、計算すべき積分を変数変換により二重指数関数的に減少する関数の全無限区間積分に変換して台形則で計算せよというものです。具体的な積分とそれに対応する変数変換は
以下のようになります。

• 区間 (−1, 1) の積分の場合
(

∫ 1

I=

−1

f (x) dx,πsinh tx = ϕDE (t) = tanh2)

• 区間 (0, ∞) の積分の場合（ただし x → ∞ で f (x) は指数関数より遅い減衰）
I=

∫ ∞

(

f (x) dx,0πsinh tx = exp
2)• 区間 (0, ∞) の積分の場合（ただし x → ∞ で f (x) は指数関数的な減衰）I=∫ ∞x = exp (t − exp(−t))

f (x) dx,0• 区間 (−∞, ∞) の積分の場合（ただし x → ±∞ で f (x) は指数関数より遅い減衰）
I=

(

∫ ∞

f (x) dx,−∞πx = sinh
sinh t2)

これら DE 公式を要約すると、減衰が二重指数関数的になるように積分形に応じて変数変換を選択することです。当然ながらもともと二重指数関数的に減衰する関数の全無限区間の積分（第 2 種変形ベッセル関数の積分表示など）はそのまま台形則で計算するのが最適になります。

5.1

DE 公式の最適性(N )高橋・森は様々な変数変換を用いたときの近似誤差 |I − Ih | の収束のオーダーを調べました [7]。

• SE 変換 ϕ(t) = tanh t の場合√(N )|I − Ih | = O(exp(−C N ))• SE+変換 ϕ(t) = tanh tm , m = 3, 5, … の場合
m

(N )

C
|I − Ih | = O(exp(− mN m+1 ))

• DE 変換 ϕ(t) = tanh( π2 sinh t) の場合(N )|I − Ih | = O(exp(−CN/ log N ))7

令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催• DE+変換 ϕ(t) = tanh( π2 sinh tm ), m = 3, 5, … の場合(N )|I − Ih | は O(exp(−CN/ log N )) より大きくなる(N )

近似誤差を離散化誤差 |I − Ih | と打切り誤差 |Ih − Ih | に分けて考えると、SE 変換
から SE+変換にすると離散化誤差は変わらずに打切り誤差が小さくなります。さらに DE 変換から DE+変換にすると離散化誤差が大きくなり全近似誤差も大きくなります。この結果から高橋・森は DE 公式が最適であると結論付けました。より数学的な最適性の証明は 1997 年に杉原によってなされました [21]。それは実軸を挟む幅 2dの帯状領域で正則な関数空間には二重指数関数より減少の速い関数は 0 以外に含まれないというものです。このことより離散化誤差が大きくなる境目がちょうど二重指数関数であることが示されました。
DE 公式の最適性は変数変換後の関数空間を適切に選んだ場合にいえます。しかし、変数変換前の関数空間をハーディ空間 H p に設定した場合には SE 変換が最適であるという事実が古くから知られています [11, 18]。これは Hp には DE 公式の多少苦手とする関数が含まれるためです。このことで DE 公式 vs SE 公式の論争が長く続き、
DE 公式が広く普及するようになったのは 1990 年後半くらいからです。

5.2

振動積分に対する DE 公式

この小節では、

∫ ∞

sin x
√
dx
0
x2 + 1のような (三角関数) × (収束の遅い関数) という形の半無限区間積分の数値計算を考えます。このような振動積分は応用上重要であるにもかかわらず、前記の DE 公式では計算困難であることが知られています [9, 15]。また古典的な補間型数値積分公式でも計算が困難です。この計算困難性の元は代数学の基本定理にあります。すなわち無限個のゼロ点が存在する振動関数は多項式では補間できないのです。この困難を克服
する 1 つの方法として、大浦・森により振動積分に対する DE 公式が提案されました[17, 19, 23, 29]。この方法はまず、目的の積分をI=∫ ∞f (x) exp(ix) dx

0

と置いて

x = M ϕ(t), ϕ(t) =で変数変換します。

I=

∫ ∞
−∞

t
1 − exp(−2t − α(1 − e−t ) − β(et − 1))f (M ϕ(t)) exp(iM ϕ(t))M ϕ0 (t) dt

ここに M , α, β は正の数とします（のちに定めます）。まだこのままでは t → +∞ での f (M ϕ(t))M ϕ0 (t) の収束が遅いので計算できません。そこで、補助的な積分

J=

∫ ∞
−∞

f (M ϕ(t)) exp(iM t)M ϕ0 (t) dt
8令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催

を導入します。そして、f (x) にある種の解析性を仮定して M → ∞ の極限を考えると |J| = O(exp(−d0 M )), (d0 > 0) のように評価されます（リーマン・ルベーグの補題の類推）。あらかじめ M を十分大きくとって |J| の値は無視できると考えて、I の代
わりに Iˆ = I − J に対して台形則を適用すると近似公式(N )Iˆh = 2iM h

N+
∑

(

)

(

)

M
M
f (M ϕ(nh)) exp iM ϕ(nh) − i ϕ̂(nh) sinϕ̂(nh) ϕ0 (nh)
2
2
n=−N−

が得られます。ここに ϕ̂(t) = ϕ(t)−t は補助的な関数です。この近似の sin(M ϕ̂(nh)/2)の部分は n → +∞ で二重指数関数的に減少するので全体として n → ±∞ で二重指数関数的に減少することになります。
(N )M = π/h となるように M と h と関連付けたとき、Iˆh の虚部だけ取り出す（f (x)を実関数として見て虚部を取り出す）とより簡単な近似∫ ∞0

f (x) sin x dx ≈ M h

N+∑f (M ϕ(nh)) sin (M ϕ(nh)) ϕ0 (nh)n=−N−で表されて、従来の変数変換型公式同様に標本点が {M ϕ(nh)}、重みが {M hϕ0 (nh)}の形になります。従来の変数変換型公式と違う点は、重みが n → +∞ で正定数に収束してゼロにはならないことです。標本点と sin の組み合わせの部分を見ると
sin(M ϕ(nh)) = (−1)n sin(M ϕ̂(nh)) となって標本点が n → +∞ で sin のゼロ点に近づくことで二重指数関数的な減少をさせていることがわかります。
(N )このときの具体的な近似誤差は h → 0 の極限で |I − Iˆh | = O(exp(−CN/ log2 N ))のように評価され、従来の DE 公式よりも少し悪い評価になります。その理由は Mが大きくなるとき f (M ϕ(t)) の極が実軸に近づき、正則な領域が狭まり離散化誤差が大きくなるためです。それを改良するため変数変換のパラメータをさらに
β = 1/4,√α = β/ 1 + M log(1 + M )/(4π) のように M と関連付けて選ぶとき、多くの場合で(N )|I − Iˆ | = O(exp(−CN/ log N )) が成り立ち、従来の DE 公式と同程度の評価が得
h

られます。

5.3

DE 公式のコンピュータへの実装

被積分関数と積分区間と許容精度を与えれば積分値と推定誤差を返すいわゆる自動
積分ルーチンを作成することを考えます [15]。DE 公式の自動積分ルーチンを作ることは比較的容易です。理由は刻み幅 h で誤差がコントロールできて誤差の漸近形もわかっているからです。例えば Ih/2 と Ih が計算できていた時の誤差を推定したい場合
は Ih の誤差はほぼ |Ih − Ih/2 | で推定できます。Ih の誤差の漸近形は A exp(−C/h) なので Ih/2 の誤差はほぼ |Ih − Ih/2 |2 /A と推定できます（A の大きさは I の複素領域積分なので |Ih | の定数倍で表されると考えます）。具体的なアルゴリズムは Ih , Ih/2 , Ih/4 ,…を誤差を推定しながら計算してゆき許容精度に達したら停止するというものです。
無限和の打ち切りについては、例えば n の絶対値が大きい時に |f (ϕ(nh))ϕ0 (nh)| の値が 2 回続けて許容精度を下回った時に打ち切るというような方法が考えられます。
ただしこれらの推定は誤る可能性があり、手動で調整しないといけない場合がありま

9令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催

す。推定の信頼度をあげる方法の 1 つとして、疑似的な積分を用いて二重の誤差推定をする方法が有効です [25]。
現在、DE 公式のルーチンは Mathematica や Maple など様々な数学ソフトウェアに組み込まれ利用されています [32]。

5.4

DE 公式の発展

前小節で紹介した積分以外にも様々な型の積分に対応する DE 公式が提案されてます [16, 20, 22, 24, 28]。
単独の DE 変換に関しては Sinc 関数近似と組み合わせることで高性能の関数近似
公式が作成できます（DE-Sinc 近似）。具体的には、区間 (−∞, ∞) で定義される関数
g(t) は()∞∑
t − nhsin πxg(t) ≈
g(nh) sinc,sinc x =hπxn=−∞で Sinc 関数近似されます。この両辺を (−∞, ∞) で積分すると全無限区間の台形則が得られることに注意します。g(t) に台形則が高精度となるときの解析性を仮定すれば、
この Sinc 関数近似も非常に高精度の近似になることがわかっています。そこで例えば DE 変換後の被積分関数 f (ϕ(t))ϕ0 (t) を g(t) と置いて Sinc 関数近似すれば高精度の被積分関数の近似ができることになり、これを応用して数値不定積分の計算ができます。この近似の多項式近似にない最大の特徴は、ある点で発散するような関数の近似ができるということです。現在、DE-Sinc 近似は不定積分の計算のほかに常微分方程式、積分方程式の近似解法などにも用いられています [26, 27, 31]。

参考文献
[1] P.J. Davis, P. Rabinowitz, Numerical integration, Blaisdell Publishing Company
1967.
[2] C. Schwartz, Numerical integration of analytic functions, J. Comput. Phys. 4(1969), 19–29.
[3] H. Takahasi, M. Mori, Error estimation in the numerical integration of analytic
functions, Rep. Compt. Centre, Univ. Tokyo 3 (1970), 41–108.
[4] 伊理正夫, 森口繁一, 高澤嘉光, ある数値積分公式について, 京都大学数理解析研
究所講究録 91 (1970), 82–118.
[5] H. Takahasi, M. Mori, Quadrature formulas obtained by variable transformation,
Numer. Math. 21 (1973), 206–219.
[6] M. Mori, On the superiority of the trapezoidal rule for the integration of periodic
analytic functions, Memoirs of Numerical Mathematics 1 (1974), 11–19.

10令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催[7] H. Takahasi, M. Mori, Double exponential formulas for numerical integration,
Publ. RIMS, Kyoto Univ. 9 (1974), 721–741.
[8] M. Mori, An IMT-type double exponential for numerical integration, Publ. RIMS,
Kyoto Univ. 14 (1978), 713–729.
[9] 戸田英雄, 小野令美, Double Exponential 変換数値積分公式の有効性を発揮させるための注意, 京都大学数理解析研究所講究録 339 (1978), 74–109.
[10] K. Murota, M. Iri, Parameter tuning and repeated application of the IMT-typetransformation in numerical quadrature, Numer. Math. 38 (1982), 347–363.
[11] K. Sikorski, F. Stenger, Optimal quadratures in H p spaces, ACM Trans. Math.
Software 10 (1984), 140–151.
[12] M. Mori, Quadrature formulas obtained by variable transformation and the DErule, J. Comput. Appl. Math. 12&13 (1985), 119–130.
[13] M. Iri, S. Moriguti, Y. Takasawa, On a certain quadrature formula, J. Comput.
Appl. Math. 17 (1987), 3–20 (English translation of [4]).
[14] M. Mori, The double exponential formulas for numerical integration over thehalf infinite interval, International Series of Numerical Mathematics 86 (1988),
367–379.
[15] 森正武, FORTRAN77 数値計算プログラミング 増補版, 岩波書店, 1988.
[16] M. Mori, Developments in the double exponential formulas for numerical integration, Proceedings of the International Congress of Mathematicians Kyoto 1990
(1991), 1585–1594.
[17] T. Ooura, M. Mori, The double exponential formula for oscillatory functions
over the half infinite interval, J. Comput. Appl. Math. 38 (1991), 353–360.
[18] F. Stenger, Numerical methods based on sinc and analytic functions, Springer
1993.
[19] M. Mori, T. Ooura, Double exponential formula for Fourier type integrals witha divergent integrand, in: R.P. Agarwal (Ed.), Contributions in Numerical Mathematics, World Scientific Series in Applicable Analysis 2 (1993), 301–308.
[20] 緒方秀教, 杉原正顯, 森正武 Cauchy の主値及び Hadamard の有限部分積分に対する DE 公式, 日本応用数理学会論文誌 3 (1993), 309–322.
[21] M. Sugihara, Optimality of the double exponential formula — functional analysis
approach —, Numer. Math. 75 (1997), 379–395.
[22] 森正武, 二重指数関数型変換のすすめ, 京都大学数理解析研究所講究録 1040
(1998), 143–153.
11令和5年度(第44回)数学入門公開講座テキスト(京都大学数理解析研究所，令和5年7月31日～8月3日開催[23] T. Ooura, M. Mori, A robust double exponential formula for Fourier type integrals, J. Comput. Appl. Math. 112 (1999), 229–241.
[24] M. Mori, M. Sugihara, The double-exponential transformation in numerical analysis, J. Comput. Appl. Math. 127 (2001), 287–296.
[25] 大浦拓哉, 二重指数関数型数値積分公式の収束判定法の改良, 日本応用数理学会
論文誌, 13 (2003), 225–230.
[26] M. Muhammad, M. Mori, Double exponential formulas for numerical indefinite
integration, J. Comput. Appl. Math. 161 (2003), 431–448.
[27] K. Tanaka, M. Sugihara, K. Murota, Numerical indefinite integration by double
exponential sinc method, Math. Comp. 74 (2004), 655–679.
[28] M. Mori, Discovery of the double exponential transformation and its developments, Publ. RIMS, Kyoto Univ. 41 (2005), 897–935.
[29] T. Ooura, A double exponential formula for the Fourier transforms, Publ. RIMS,
Kyoto Univ. 41 (2005), 971–978.
[30] T. Ooura, An IMT-type quadrature formula with the same asymptotic performance as the DE formula, J. Comput. Appl. Math. 213 (2008), 232–239.
[31] 田中健一郎, 岡山友昭, 変数変換型数値計算法, 岩波書店, 2023.
[32] Wolfram 言語ドキュメント―NIntegrate 積分ストラテジー,
https://reference.wolfram.com/language/tutorial/NIntegrateIntegrationStrategies.html

12