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

振動積分の例
振動積分の例: The SIAM 100-Dollar, 100-Digit Challenge (2002)∫ 1

Problem 1:

What is lim

ε→0 ε

x−1 cos(x−1 log x)dx ?

30x-1cos(x-1logx)

20100-10-20

-30
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1x

・従来の方法では計算が困難
・理工学の分野で重要な計算問題
1

従来の振動積分の計算法
従来の方法:

I ==

∫ ∞
0∞
∑f (x) sin x dx
∫ (n+1)πSn,n=0Sn =nπf (x) sin x dx

ステップ 1 S0, S1, S2, · · · , SN の値を計算 (積分計算)∑∞
ステップ 2 n=0 Sn の値を S0, S1, S2, · · · , SN から推測する (補外)補外法の例: オイラー変換、チェビシェフ加速、エイトケンの ∆2 加速など
この方法の問題点:
・たくさんの関数計算が必要
・誤差の見積もりが困難
2

振動積分に対する DE 公式の歴史• 1974: DE 公式の提案 (高橋・森)
この DE 公式は、広義積分が計算できる万能型数値積分公式だが、振動積分は計算できない。
• 1991: 振動型 DE 公式の提案 (大浦・森)
1989 年に、 私が大学生だったころ、この方法を発見し森正武氏に手紙を書くことから始まる。
• 1999: 振動型 DE 公式の改良 (大浦・森)• 2005: 振動型 DE 公式の一般化 (大浦)3振動型 DE 公式の導出F =∫ ∞f (x) exp(iωx) dx,

ω>0w 0w DE 変換: x = M φ(t),wwtφ(t) = 1−exp(−2t−α(1−e−t )−β(et −1)) ,∫ ∞F =f (M φ(t)) exp (iωM φ(t)) M φ0(t) dt−∞

M, α, β > 0

擬似的な積分:E=∫ ∞−∞()f (M φ(t)) exp iωM φ(t) − iω0M φ̂(t) M φ0(t) dt,

ここで φ̂(t) = φ(t) − t
ベーグの補題の類推）⇒|E| ≈ 0 for large M （リーマン・ルF の代わりに F̃ = F − E を計算する。
4

F , E, F̃ の被積分関数
Integrand of F

420-2-4-4-3-2-1012t

Integrand of F-E

Integrand of E

420-2-4-4-3-2-10t1286420-2-4-6-8-4-3-2-101t52

F̃ の収束性F̃ =∫ ∞()

if (M φ(t)) exp iωM φ(t) − ω0M φ̂(t)2−∞()
1· 2iM sin
ω0M φ̂(t) φ0(t) dt2被積分関数の収束:
・φ0(t) = o(exp(−αe|t|)) (t → −∞)・φ̂(t) = φ(t) − t = o(exp(−βe|t|)) (t → +∞)()1⇒ sin 2 ω0M φ̂(t) φ0(t) → 0 Double Exponentially as t → ±∞6

振動型 DE 公式
N+∑2πi(N )F̃h=fω0 n=−N(−· φ0(nh) exp)()

ππφ(nh) sinφ̂(nh)ω0h2h()πiω
πiφ(nh) −φ̂(nh)ω0 h2hwheretφ(t) =,−tt
1 − exp(−2t − α(1 − e ) − β(e − 1))√α = β/ 1 + log(1 + π/(ω0h))/(4ω0h), β = 0.25.
ω の範囲と誤差:(N )|F − F̃h0 < ω < 2ω0,| < c00e−c0/h + c01e−c1ω/h + c02e−c2(2ω0−ω)/h< C 0(ω) exp(−C(ω)N/ log N )7

振動型 DE 変換の変遷• 1989: 最初の提案（未発表）√t + t2 + A exp(−B cosh t)
φ(t) =2A, B > 01

• 1989,1991: 最初の論文での提案

0.5tφ(t) =1 − exp(−6 sinh t)

0-101• 1999,2005: 改良型の変換

tφ(t) =
1 − exp(−2t − α(1 − e−t) − β(et − 1))h )1/2 ), β = O(1)α = O(( − logh8計算例
例 1: F =∫ ∞0log x√ exp(iωx) dxx

パラメータの選択: ω0 = 1, h = 0.075, N− = 94, N+ = 69100
10-2Error

10-4
10-610-810-10
10-12
10-14
10-1600.5

11.52

ω(N )

誤差 |F − F̃h|9計算例 (発散振動積分)
例 2: F (ω) =∫ ∞0log x exp(iωx) dx，F = lim F (ω + iε) で定義
ε→0

パラメータの選択: ω0 = 1, h = 0.075, N− = 87, N+ = 65
10010-2Error

10-410-610-810-10
10-1210-14
00.5

11.52

ω

(N )

誤差 |F − F̃h|10

高精度計算例 (SIAM-100digit)∫ 1
0x−1 cos(x−1 log x)dx =∫ ∞0cos u
du−1u(1 + (W (u)) )(W :ランベルトの W 関数)

パラメータ：ω0/2 = ω = 1
標本点数刻み幅 h
誤差381/8
2.2 · 10−11901/16
1.5 · 10−20207
1/32
1.4 · 10−38

471
1/64
7.5 · 10−73

1057
1/128
1.8 · 10−14311

DE 公式の自動積分ルーチン（概略）∫ 1

I =f (x) dx

−1∞
∑(())

π cosh(nh)2 ()π2cosh 2 sinh(nh)))((π cosh(nh)π2|
無限和の打ち切り：指定の精度 ε よりも |f tanh 2 sinh(nh)cosh2 ( π2 sinh(nh))πf tanh
Ih = hsinh(nh)2n=−∞が小さくなったところで和を打ち切る。

h の決定アルゴリズム：|I − Ih| ≤ C 0||f || exp(−C/h), C, C 0 > 0 から刻み幅 h を半分にすると、有効桁数が倍になることを利用してIh, Ih/2, Ih/4, Ih/8, · · ·を誤差推定を行いながら計算し，誤差推定が ε を下回ったとことで終了する。
12離散化誤差の推定（粗い判定）
仮定：h = h0, h0/2, h0/22, · · · , h0/2M に対して

α|I| exp(−C/h) ≤ |I − Ih| ≤ β|I| exp(−C/h)となる α, β > 0 が存在する。

Eh = |I − Ih|/|I| ≈ |Ih/2 − Ih|/|Ih/2|
Eh/2 ≤ β exp(−2C/h) ≤β 2Eh2α

収束判定：ε を許容相対誤差とするとき

α √ε|Ih/2 − Ih|/|Ih/2| < √βを満たしたとき計算を終了する。
√
ここで α/ β ：安全因子
13離散化誤差の推定の問題点
安全因子の定め方：いくつかの積分を用意して経験的に定める。
問題点：ある有理関数の積分の正確な誤差評価

I − Ih ≈ < [A exp(−B/h)]= |A| cos(argA − =B/h) exp(−<B/h)h が cos のゼロ点となる値になった場合、α は極端に小さくなる。
したがって安全因子が小さくなり収束判定が厳しくなるので計算量が増える。

14離散化誤差 |I − Ih| の図
010-21010-4|Ih - I|

10-610-810-1010-1210-1410-16

12345671/h∫ 1

I=−1(1.2 + x)−1 dx の離散化誤差15重み付き積分の同時計算 (誤差推定の改良)∞∑(())∞∑(())

πf tanhsinh(nh)Ih = h2n=−∞π cosh(nh)2 ()π2cosh 2 sinh(nh)と同時に

πJh = hf tanhsinh(nh)2n=−∞(2)π2cosh 2 sinh(nh)を計算する。これは，重みつき積分
∫ 1

J=−1

√

4

π 2 + 4(tanh−1 x)2

f (x) dxを計算することに相当Jh の離散化誤差：|J − Jh| ≤ C 00||f || exp(−C/h),

C 00 > 016

離散化誤差の推定 (改良)仮定 (改良)：h = h0, h0/2, h0/22, · · · , h0/2M に対して

α0|I| exp(−C/h) ≤ max(|I − Ih|, |J − Jh|),|I − Ih| ≤ β|I| exp(−C/h)となる α0, β > 0 が存在する。

Eh0 = max(|I − Ih|, |J − Jh|)/|I|≈ max(|Ih/2 − Ih|, |Jh/2 − Jh|)/|Ih/2|
β 02
Eh/2 ≤ β exp(−2C/h) ≤ 02 Ehα
収束判定 (改良)：以下の式を満たしたとき計算を終了する。
α0 √εmax(|Ih/2 − Ih|, |Jh/2 − Jh|)/|Ih/2| < √β17max(|Ih/2 − Ih|, |Jh/2 − Jh|) の図
010-2max(|Ih - I|, |Jh - J|)

1010-410-6
10-8
10-10
10-12
10-1410-16

12345671/h∫ 1I=−1(1.2 + x)−1 dx の離散化誤差18指数関数計算をほぼ含まない DE 公式— DE 変換の計算における指数関数の排除 —アイデア：
1. 刻み幅を h = log 2/K, (K は整数) に選ぶ．
2. 一番内側のループを刻み幅 Kh ステップで計算する．
⇒ 加減算と乗算だけで重みと標本点が計算可能∵ an+1 = exp(exp((n + 1)Kh)) = exp(2 exp(nKh)) = a2n
ただし、初期計算に log 2 の計算と少しの指数関数と平方根の計算が必要になる。

19

初等関数計算の高速化
よく知られた初等関数の算法：
・CORDIC 法：N 桁計算に O(N ) 回の加算とビットシフトが必要
・テイラー展開による方法：N 桁計算に O(N ) 回の乗算が必要√
・加法定理で区間を狭めてテイラー展開：N 桁計算に O( N ) 回の乗算・ガウスの算術幾何平均に基づく方法：N 桁計算に O(log N ) 回の乗算
高精度高速算法の例：
・CORDIC 法で区間を狭めてテイラー展開：√
N 桁計算に O( N ) 回の乗算，加算とビットシフトが必要√
計算量の N にかかる比例係数が小さい
巨大なテーブルを用いることで徹底した最適化が可能
20

より高性能の DE 変換の例)

(πsinh t : 高橋・森 (1974)
ϕDE(t) = tanh2よりも性能の良い変換∫t1
ϕDEa1(t) =exp(−k cosh s) ds : 杉原 (1997)
K0(k) 0(())π4sinh t: 大浦 (2016)
ϕDEa2(t) = −1 + tan−1 exp
π2
特に ϕDEa2(t) は ϕDE(t) よりも 2 割から 3 割程度計算量が少なくなる。

21

DE 変換のさらなる発展・コーシの主値積分など様々な積分に関する DE 公式が提案されている。
・DE 変換と Sinc 関数近似の組み合わせで関数近似ができる (DE-Sinc 近似)。
・DE-Sinc 近似を用いて数値不定積分の計算や常微分方程式の数値解法などが提案されている。
これらの DE 変換に関する算法の特徴：
・高速かつ高精度である。
・端点の特異性に強い。
22