数値解析と複素関数論 (7)
数値積分：DE 公式
緒方 秀教
電気通信大学

大学院情報理工学研究科 情報・ネットワーク工学専攻

2021 年 4 月 13 日（火）緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

前回までの復習：全無限区間積分に対する台形則
台形則は解析関数の全無限区間積分に強い．
Z ∞
N+∞
XX
f (kh)f (kh) ≃ hf (x)dx ≃ h
−∞

高精度.

k=−N−k=−∞f (z) は帯状閉領域
Dd ≡ { z ∈ C | | Im z| ≤ d }(d >0)を含む領域で解析的である．

N(f , d) :=I|f (z)||dz| < ∞,∂Ddlimx→±∞

Z d|f (x + iy )|dy = 0.

−d

このとき，台形則は指数関数的収束する：

Z ∞−∞f (x)dx − h
2πd.
f (kh) . N(f , d)exp −hk=−∞∞X緒方秀教数値解析と複素関数論 (7) 数値積分：DE 公式

前回までの復習：変数変換を用いた数値積分
台形則は解析関数の全無限区間積分に強い．

→ この事実を任意区間の積分に応用したい．
変数変換 x = ψ(u) の利用．

Z baf (x)dx =Z ∞−∞f (ψ(u))ψ ′ (u)du ≃ hN+X

f (ψ(kh))ψ ′ (kh).

k=−N−

変数変換の例：SE 変換（一重指数関数的変数変換）

1/2,
cosh2 (u/2)′ψSE(u) ≃ 2 exp(−|u|) ( u → ±∞, 一重指数関数的減衰 ).
uψSE (u) = tanh ,2′ψSE(u) =少ない標本点数 N+ + N− + 1 で台形則の無限和を打ち切れる．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

今回の内容もっと効率の良い変数変換はないか？

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

今回の内容もっと効率の良い変数変換はないか？

→ DE 公式（二重指数関数型公式，double exponential formula）

H. Takahasi & M. Mori: Double exponential formulas for numerical integration,
Publ. RIMS Kyoto Univ., 9 (1974) 721–741.

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

今回の内容もっと効率の良い変数変換はないか？

→ DE 公式（二重指数関数型公式，double exponential formula）

H. Takahasi & M. Mori: Double exponential formulas for numerical integration,
Publ. RIMS Kyoto Univ., 9 (1974) 721–741.

今回の内容
1DE 公式（公式，数値例）

2

DE 公式の最良性

3

DE 公式の使用上の注意

4

無限区間積分に対する DE 公式

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

1/4 DE 公式

DE 公式（二重指数関数型公式，double exponential formula）
Z 1−1f (x)dx ≃ h

ψDE (u) = tanh(c sinh u)

N+X′f (ψDE (kh))ψDE(kh),k=−N−( c > 0, DE 変換（二重指数関数型変換）)

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

1/4 DE 公式
DE 公式（二重指数関数型公式，double exponential formula）
Z 1−1f (x)dx ≃ h−1f (x)dx ≃ h≃h′f (ψDE (kh))ψDE(kh),k=−N−( c > 0, DE 変換（二重指数関数型変換）)

ψDE (u) = tanh(c sinh u)

Z 1N+X∞X′f (ψDE (kh))ψDE(kh) もともとは無限和

N+X′f (ψDE (kh))ψDE(kh).

k=−∞k=−N−′|f (ψDE (kh))ψDE(kh)| < （ǫ 非常に小さい数）となったところで無限和を打ち切る．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

1/4 DE 公式
DE 公式（二重指数関数型公式，double exponential formula）
Z 1−1f (x)dx ≃ h

ψDE (u) = tanh(c sinh u)

N+X′f (ψDE (kh))ψDE(kh),k=−N−( c > 0, DE 変換（二重指数関数型変換）)
1.5c cosh ucosh2 (c sinh u)≃ 2c exp(|u| − c exp |u|)′ψDE(u) =( u → ±∞ ).

10.5
0-0.5
psi(u)psi’(u)-1-4′
ψDE(u) は二重指数関数的減衰する．

-3-2-10u1

234( c = π/2 )→ もっと少ない標本点数 N+ + N− + 1 で台形則の無限和を打ち切れる．
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

1/4 DE 公式：数値例
Z 1πdx= = 1.57079 63267 94896 61923 13216 91640 . . .
22−1 1 + x

多倍長（10 進 200 桁）演算．|f (kh)| < 10−200 で無限和打ち切り．
11DESE
1x10
relative error

relative error

1×10DESE−50

1×10−100
1×10−150
1×10−200-50

1x10-100
1x10-150

0204060 80 100 120 140
1/h

1x10-200020000

40000
N

相対誤差のグラフ
横軸：1/h横軸：標本点数

DE 変換 : ψDE (u) = tanh(c sinh u) ( c = π/2 ),SE 変換 : ψSE (u) = tanh(u/2).
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

60000

1/4 DE 公式：数値例
Z 1πdx= = 1.57079 63267 94896 61923 13216 91640 . . .
22−1 1 + x多倍長（10 進 200 桁）演算．|f (kh)| < 10−200 で無限和打ち切り．
11DESE
1x10
relative error

relative error

1×10DESE−50

1×10−100
1×10−150
1×10−200-50

1x10-100
1x10-150

0204060 80 100 120 140
1/h

1x10-200

02000040000
N

相対誤差のグラフ
横軸：1/h横軸：標本点数

DE 公式は，SE 公式に比べてはるかに少ない標本点数で積分値を計算できる．
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

60000

1/4 DE 公式

DE 公式は端点特異性をもつ積分に有効である．
Z 1
Z 1
f (x)dx =f0 (x)(1 + x)α−1 (1 − x)β−1−1( α, β > 0 ).

−1

緒方

秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

1/4 DE 公式

DE 公式は端点特異性をもつ積分に有効である．
Z 1Z 1f (x)dx =f0 (x)(1 + x)α−1 (1 − x)β−1−1( α, β > 0 ).

−1f (x) = f0 (x)(1 − x 2 )−1/2 の場合
Z 1−1f (x)dx =

Z ∞

−∞

c cosh u
f0 (tanh(c sinh u)) (1 − tanh2 (c sinh u))−1/2
du{z} cosh2 (c sinh u)|cosh(c sinh u)

Z ∞c cosh udu.
=f0 (tanh(c sinh u))cosh(csinh u)−∞

f (x) が端点 x = ±1 にベキ的特異性をもっても，DE 変換変換後の被積分関数は指数関数的減衰する．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

1/4 DE 公式

DE 公式は端点特異性をもつ積分に有効である．
Z 1
Z 1f (x)dx =f0 (x)(1 + x)α−1 (1 − x)β−1−1( α, β > 0 ).

−1

一般の場合

Z 1f (x)dx =−1=Z 1

f0 (x)(1 + x)α−1 (1 − x)β−1 dx−1Z ∞f0 (ψDE (u))−∞

c cosh u exp(c(α − β) sinh u)
coshα+β (c sinh u).

f (x) が端点 x = ±1 にベキ的特異性をもっても，DE 変換変換後の被積分関数は指数関数的減衰する．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

1/4 DE 公式：数値例
Z 1−1√dx= π = 3.14159 26535 89793 23846 26433 83280 . . .
1 − x2

多倍長（10 進 200 桁）演算．|f (kh)| < 10−200 で無限和打ち切り．
11DESE
1×10relative error

relative error

1×10DESE−50

1×10−100
1×10−150
1×10−200−50

1×10−100
1×10−150
1×10−200

020401/h

6080020000

40000
N

相対誤差のグラフ
横軸：1/h横軸：標本点数

DE 変換 : ψDE (u) = tanh(c sinh u) ( c = π/2 ),SE 変換 : ψSE (u) = tanh(u/2).
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

60000

1/4 DE 公式：数値例
Z 1−1√dx= π = 3.14159 26535 89793 23846 26433 83280 . . .
1 − x2

多倍長（10 進 200 桁）演算．|f (kh)| < 10−200 で無限和打ち切り．
11DESE
1×10relative error

relative error

1×10DESE−50

1×10−100
1×10−150
1×10−200−50

1×10−100
1×10−150
1×10−200020401/h

6080020000

40000
N

相対誤差のグラフ
横軸：1/h横軸：標本点数

DE 公式は，SE 公式に比べてはるかに少ない標本点数で積分値を計算できる．
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

60000

2/4 DE 公式の最良性

DE 変換は最良な変数変換である．
ψ ′ (u) がこれ以上速く減衰するような変数変換（三重指数関数型変換，…）を使っても，効率は上がらない．
理論的証明．

M. Sugihara: Optimality of the double exponential formula —functional analyris approach —, Numer. Math. 75 (1997) 379–395.

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

2/4 DE 公式の最良性

DE 変換は最良な変数変換である．
「三重指数関数型変換」をやってみても精度は上がらない．
1DE 変換およびDETETE（三重指数関数型）変換-50
relative error

1x10x = tanh(c sinh(sinh u))を用いて積分

1x10-100
1x10-150
1x10-200

Z 1dx√
1 − x2−10200

400
N600800

縦軸： 相対誤差
横軸： 標本点数を計算した結果．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

2/4 DE 公式の最良性
DE 変換が最良であることの直感的理解．
ψ ′ (u) が急減衰するほど，標本点ψ(kh) は積分端点 x = ±1 の近くに密集する．
あまり ψ ′ (u) を急減衰させると，積分区間 (−1, 1) の中央あたりの標本点分布がスカスカになる．

→ 被積分関数 f (x) の情報がうまく拾えない．
何処かに ψ ′ (u) の急減衰の速さが最適なところがある．

. . . それが DE 変換．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

3/4 DE 公式の使用上の注意
端点特異性を持つ積分を DE 公式で計算する場合，桁落ち対策をしなければならない．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

3/4 DE 公式の使用上の注意
端点特異性を持つ積分を DE 公式で計算する場合，桁落ち対策をしなければならない．

Z 1−1f0 (x)(1 − x 2 )−1/2 dx =

Z ∞

−∞

≃h′f0 (ψDE (u))(1 − ψDE (u)2 )−1/2 ψDE(u)duN+X

k=−N−′f0 (ψDE (kh))(1 − ψDE (kh)2 )−1/2 ψDE(kh),ψDE (kh) = tanh(c sinh kh) → ±1 ( k → ±∞, 急接近 ).
|k| が大きいとき，
1 − ψDE (kh)2 = 1 − (1 に非常に近い数)= (有効桁数の小さい，非常に小さい数).

(1 − ψDE (kh)2 )−1/2 は有効桁数の小さい非常に大きい数となり，数値積分の精度が落ちる．
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

3/4 DE 公式の使用上の注意
Z 1=f0 (x)(1 − x 2 )−1/2 dx

−1Z ∞

−∞

=Z ∞

−∞

Z 1−1f0 (ψDE (u))(1 − tanh2 (c sinh u))−1/2f0 (ψDE (u))

c cosh u
ducosh2 (c sinh u)c cosh udu,cosh(c sinh u)

f0 (x)(1 − x 2 )−1/2 du ≃ hN+X

k=−N−f0 (ψDE (kh))

c cosh(kh). (1)cosh(c sinh(kh))(1) 右辺の計算では桁落ちが生じない．
あらかじめ積分を (1) 右辺の形に手計算で直してから，数値計算をする．
コンピュータに “1 − (1 に非常に近い数)”の計算をさせない．
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

3/4 DE 公式の使用上の注意

一般の場合

Z 1−1=f0 (x)(1 + x)α−1 (1 − x)β−1 dx

Z ∞f0 (ψDE (u))

c cosh u exp(c(α − β) sinh u)
coshα+β (c sinh u)

−∞

≃hN+X

f0 (ψDE (kh))

duc cosh kh exp(c(α − β) sinh kh)k=−N−coshα+β (c sinh kh)

コンピュータプログラムには最右辺の形で数値積分公式を書く．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

（補足事項）DE 変換の最適パラメータ
DE 変換 ψDE (u) = tanh(c sinh u) の係数 c ．
c=πが最適である．
2（理由）全無限区間積分に対する台形則．
Z ∞∞Xg (kh).
g (u)du ≃ h−∞

k=−∞


2πdg (w ) が | Im w | ≤ d で正則，etc. ならば，誤差 = O exp −.
h
帯幅 2d が大きいほど（g (w ) が正則である範囲が広いほど）精度が良い．
′g (w ) = f (ψDE (w ))ψDE(w ) が正則である帯状領域 | Im w | ≤ d ができるだけ広くなるように c をとると，c = π/2 が最適であることがわかる．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

（補足事項）DE 公式の理論誤差評価
Z 1f (x)dx − h

Z 1f (x)dx − h−1≤−1N0X′f (ψDE (kh))ψDE(kh)∞X′f (ψDE (kh))ψDE(kh) + h. . . DE 公式の誤差

k=−N0k=−∞X′f (ψDE (kh))ψDE(kh)|k|>N0=: 離散化誤差 + 打切り誤差.



2πd
離散化誤差：O exp −h.

打切り誤差：O[exp(−α exp Nh)](α > 0 : const.)

離散化誤差 ≈ 打切り誤差，となるように，標本点数 N(= 2N0 + 1)
に対し刻み幅 h を取れば，

N.
DE 公式の誤差 = O exp −const. ×log N
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

4/4 無限区間積分に対する DE 公式
Z ∞0f (x)dx =Z ∞−∞f (ψ(u))ψ ′ (u)du ≃ h緒方秀教

N+Xf (ψ(kh))ψ ′ (kh).

k=−N−

数値解析と複素関数論 (7) 数値積分：DE 公式

4/4 無限区間積分に対する DE 公式
Z ∞0f (x)dx =Z ∞−∞f (ψ(u))ψ ′ (u)du ≃ hf (x) がベキ的減衰：f (x) = O(|x|−α )N+Xf (ψ(kh))ψ ′ (kh).

k=−N−( x → +∞, α > 1 ).

ψ(u) = exp(2 sinh u).
f (x) が指数関数的減衰：f (x) = O[exp(−x)]( x → +∞ ).

ψ(u) = exp(u − exp(−u)).
∗ 前項の変換では被積分関数が三重指数関数的減衰になってしまう．

f (x) が Gaussian で減衰：f (x) = O[exp(−x 2 )] ( x → +∞ ).
u− exp(−u) .
ψ(u) = exp2緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式4/4 無限区間積分に対する DE 公式
減衰の遅い振動関数の無限区間積分（Fourier 変換）に対しては（上記の）DE 公式は有効でない．
Z ∞f0 (x) sin(ωx + β)dx, f0 (x) = O(|x|−α ) ( x → +∞ ).
0緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

4/4 無限区間積分に対する DE 公式
減衰の遅い振動関数の無限区間積分（Fourier 変換）に対しては（上記の）DE 公式は有効でない．
Z ∞f0 (x) sin(ωx + β)dx, f0 (x) = O(|x|−α ) ( x → +∞ ).
0

振動積分に特化した DE 公式．

T. Ooura & M. Mori: The double exponential formula for
oscillatory functions over the half infinite interval,
J. Comput. Appl. Math., 38 (1991) 353–360.
T. Ooura & M. Mori: A robust double exponential formula for
Fourier type integrals, J. Comput. Appl. Math., 112 (1999)
229–241.
佐藤超函数論的手法．

H. Ogata: Numerical calculation of Fourier transforms basedon hyperfunction theory, J. Comput. Appl. Math., 378 (2020)112921.
緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式

まとめDE 公式
DE 変換：被積分関数が二重指数関数的減衰するような変数変換．
端点特異性を持つ積分に対しても有効である．
最適な変数変換．
DE 公式は，端点特異性を持つ積分を計算する場合，桁落ち対策をしなければならない．
無限区間積分に対する DE 公式．

緒方秀教

数値解析と複素関数論 (7) 数値積分：DE 公式