using DoubleExponentialFormulas
using QuadGK
using Printf
using Plots
using LaTeXStrings
using SpecialFunctions
基本的に数値積分にかかる時間は、被積分関数の評価にかかる時間が大半を占めます。なので、実行時間を計っても被積分関数を評価するばかりで、数値積分関数自体の評価のためにはよい指標とは言えません。かわりに被積分関数 $f(x)$ を評価した回数を指標として用いることが一般的なようです。目標とする精度の積分値をより少ない回数で出せるほうがよい数値積分法といえます。
注: ただし、目標とする精度ぴったりを目指して計算することはほとんどの場合困難です。通常は、適当に反復をしながら目標精度を上回った時点で打ち切る、という方法をとるので反復のパラメータを微調整することで性能が劇的に改善することがあります。
注: ひとまず、以下では quadde
, quadgk
ともにデフォルトの精度での計算を行っていきます。これらは一致しており、およそ7~8桁の数字が厳密解と一致する程度を目標としています。
以下の @printeval
マクロは数値積分を実行し、被積分関数の評価回数、積分値、推定誤差を画面に表示します。 EvalCounter
は @printeval
マクロの内部で使われる型です。
mutable struct EvalCounter <: Function
n::Int
f::Function
EvalCounter(f::Function) = new(0, f)
end
function (ec::EvalCounter)(x)
ec.n += 1
ec.f(x)
end
macro printeval(expr)
s = string(expr)
q = expr.args[1]
f = expr.args[2]
rest = expr.args[3:end]
quadexpr = Expr(:call, esc(q), :ec)
for part in rest
if part isa Expr || part isa Symbol
push!(quadexpr.args, esc(part))
else
push!(quadexpr.args, part)
end
end
blk = Expr(:block)
push!(blk.args, :(ec = EvalCounter($(esc(f)))))
push!(blk.args, Expr(:(=), :ret, quadexpr))
push!(blk.args, :(println($(s) * ": Evaluation=", @sprintf("%-5d", ec.n),
" I=", @sprintf("%-1.15e", ret[1]),
" E=", @sprintf("%-1.4e", ret[2]))))
push!(blk.args, :ret)
return blk
end
@printeval (macro with 1 method)
これらは明らかに不定積分が可能なので数値積分を使う必要はありませんが、ひとまず例として挙げます。 quadgk
で使われているガウスークロンロッド積分法は低次の多項式に関しては、原理的にほぼ厳密解を返します(うろおぼえ)。 quadgk
は驚異的に少ない評価回数で解を返しています。 quadde
も二桁以内の回数に収めているので十分だと思います。
注: これら低次の多項式に限らず、以下で紹介する積分は複雑に見えても厳密解を得られるものが多いです。これは性能評価において、まず数値積分が正しく動いているかどうかを確かめる必要があるためです。しかし、現実に数値積分が必要な場合においては厳密解を得ることは不可能で、被積分関数 $f(x)$ はより複雑であると考えられます。
let f(x) = x
@printeval quadde(f, 0.0, 1.0)
@printeval quadgk(f, 0.0, 1.0)
plot(f, label=L"x", xlim=[0.0, 1.0],
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, 1.0): Evaluation=25 I=5.000000000000183e-01 E=1.1287e-10 quadgk(f, 0.0, 1.0): Evaluation=15 I=5.000000000000000e-01 E=0.0000e+00
let f(x) = x^3 + x^2 + x + 1
@printeval quadde(f, -1.0, 1.0)
@printeval quadgk(f, -1.0, 1.0)
plot(f, label=L"x^3 + x^2 + x + 1", xlim=[-1.0, 1.0],
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, -1.0, 1.0): Evaluation=51 I=2.666666666666667e+00 E=5.9246e-16 quadgk(f, -1.0, 1.0): Evaluation=15 I=2.666666666666667e+00 E=0.0000e+00
let f(x) = exp(-x^2)
@printeval quadde(f, -Inf, Inf)
@printeval quadgk(f, -Inf, Inf)
plot(f, label=L"\exp{\left(-x^2\right)}", xlim=[-5.0, 5.0],
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, -Inf, Inf): Evaluation=227 I=1.772453850904036e+00 E=2.5628e-12 quadgk(f, -Inf, Inf): Evaluation=225 I=1.772453850905514e+00 E=6.4296e-09
let f(x) = x^3*exp(-x^2)
@printeval quadde(f, 0.0, Inf)
@printeval quadgk(f, 0.0, Inf)
plot(f, label=L"x^3 \exp{\left(-x^2\right)}", xlim=[0.0, 5.0],
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, Inf): Evaluation=145 I=4.999999999999333e-01 E=2.3797e-12 quadgk(f, 0.0, Inf): Evaluation=165 I=5.000000000000004e-01 E=2.7808e-09
理由はわかりませんが、被積分関数に平方根の計算 sqrt
が含まれると数値積分の性能が悪化しやすくなるような気がします。(有効桁数の問題?)しかし、以下二つの簡単な例ではなぜか quadde
が非常にいい結果を出しています。
let f(x) = sqrt(x)
@printeval quadde(f, 0.0, 1.0)
@printeval quadgk(f, 0.0, 1.0)
plot(f, label=L"\sqrt{x}", xlim=[0.0, 1.0],
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, 1.0): Evaluation=25 I=6.666666666668387e-01 E=2.6111e-09 quadgk(f, 0.0, 1.0): Evaluation=315 I=6.666666670773991e-01 E=7.1095e-09
let f(x) = 1/sqrt(x)
@printeval quadde(f, 0.0, 1.0)
@printeval quadgk(f, 0.0, 1.0)
plot(f, label=L"\frac{1}{\sqrt{x}}", xlim=[0.0, 1.0],
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, 1.0): Evaluation=25 I=1.999999982931068e+00 E=9.0356e-12 quadgk(f, 0.0, 1.0): Evaluation=1305 I=1.999999984598392e+00 E=2.3763e-08
一般的に言うと、振動する関数は数値積分が少し困難です。特に周期が極端に短い場合は非常に困難になり、積分区間を分割するなどの対策が必要になります。しかし、以下の簡単な例では、うまくいっているようですね。 quadgk
の評価回数はここでも驚異的に少ないです。
let f(x) = sin(x)
@printeval quadde(f, 0.0, π)
@printeval quadgk(f, 0.0, π)
plot(f, label=L"\sin{x}",
xlim=[0.0, π], xticks=(0.0:π/2:π, ["0", "\\pi/2", "\\pi"]),
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, π): Evaluation=51 I=2.000000000000000e+00 E=1.1712e-14 quadgk(f, 0.0, π): Evaluation=15 I=2.000000000000000e+00 E=1.7906e-12
let f(x) = cos(x)
@printeval quadde(f, 0.0, π/2)
@printeval quadgk(f, 0.0, π/2)
plot(f, label=L"\cos{x}",
xlim=[0.0, π/2], xticks=([0.0, π/4, π/2], ["0", "\\pi/4", "\\pi/2"]),
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, π / 2): Evaluation=51 I=1.000000000000000e+00 E=2.2393e-16 quadgk(f, 0.0, π / 2): Evaluation=15 I=1.000000000000000e+00 E=2.2204e-16
let f(r) = exp(-(r^2 - 1))
y = 1
g(r) = 2*f(r)*r/sqrt(r^2 - y^2)
@printeval quadde(g, y, Inf)
@printeval quadgk(g, y, Inf)
plot(g, label=L"\frac{2 f(r) r}{\sqrt{r^2 - y^2}}",
xlim=[y, 2.0], fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(g, y, Inf): Evaluation=89 I=1.772453810111211e+00 E=5.3962e-11 quadgk(g, y, Inf): Evaluation=1395 I=1.772453829887660e+00 E=4.3196e-09
もとの式は $r = y$ となるすべての点で発散しますが、 $r = s^2 + |y|$ (y は実数)とおいて変数変換すると $r = y = 0$ の一点のみで発散する式に改良できます。
$$ F(y) = 2 \int_{|y|}^{\infty} \frac{f(r) r}{\sqrt{r^2 - y^2}} dr = 2 \int_{0}^{\infty} \frac{f(s^2 + |y|) (s^2 + |y|)}{\sqrt{\left( s^2 + |y| \right)^2 - y^2 }} \frac{dr}{ds} ds = 4 \int_{0}^{\infty} \frac{f(s^2 + |y|) (s^2 + |y|)}{\sqrt{s^2 + 2|y|}} ds = 4 \int_{0}^{\infty} \frac{f(r) r}{\sqrt{r + |y|}} ds $$被積分の変形によって quadgk
の評価回数が 1/10 に改善されました。(quadde
はなぜか変換後の関数のほうが苦手のようです)
let f(r) = exp(-(r^2 - 1))
y = 1
function g(s)
absy = abs(y)
r = s^2 + absy
4*f(r)*r/sqrt(r + absy)
end
@printeval quadde(g, 0.0, Inf)
@printeval quadgk(g, 0.0, Inf)
plot(g, label=L"\frac{4 f(s^2 + |y|) (s^2 + |y|)}{\sqrt{s^2 + 2|y|}}",
xlim=[0.0, 2.0], fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(g, 0.0, Inf): Evaluation=161 I=1.772453850904258e+00 E=9.1114e-13 quadgk(g, 0.0, Inf): Evaluation=135 I=1.772453850905516e+00 E=4.0140e-09
ほとんどの数値積分アルゴリズムは連続な関数 $f(x)$ の積分を離散値の重み付き和に近似させて計算します。
$$ \int_{a}^{b} f(x) dx \approx \Sigma^{k} f(x_k) w_k $$ここで $x_k$ は離散値であり、非適応型の積分法では $x_k$ を一定の間隔でとります。結果的に、最も変化の激しい(=スパイク)部分にあわせて全体を細かい間隔で評価するために、このような被積分関数に対しては $k$ の最大値が大きくなる場合があります。対して、適応型の積分法では $f(x)$ の変化の激しさによって $x_k$ の間隔を調整するためこの問題は幾分緩和されます。しかし、スパイクに含まれる $x_k$ を「発見」できるか、という困難は常に存在しています。このスパイクを発見できるか、という問題は適応型・非適応型問わずに存在する問題で、ユーザーが結果を見つつ対応しなければなりません。
$$ \int_{0}^{1} \cosh^{-2}{(10x-2)} + \cosh^{-4}{(100x-40)} + \cosh^{-6}{(1000x-600)} dx = 0.2108027355005492773756\dots $$上の積分は x = 0.4
および x = 0.6
付近にスパイクを含みます。適応型の quadgk
は評価回数を抑えていますが 0.01 程度の誤差があります、どうやら一つのスパイクを見落としているようです。対して、非適応型の quadde
は9桁を合わせていますが、評価回数は5桁に及んでしまっています。
let f(x) = 1/cosh((10x-2))^2+1/cosh((100x-40))^4+1/cosh((1000x-600))^6
@printeval quadde(f, 0.0, 1.0)
@printeval quadgk(f, 0.0, 1.0)
x = 0.0:0.001:1.0
plot(x, f.(x), label="", xlim=[0.0, 1.0], xticks=0.0:0.2:1.0,
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, 1.0): Evaluation=12893 I=2.108027351307217e-01 E=3.2400e-10 quadgk(f, 0.0, 1.0): Evaluation=315 I=2.097360711575962e-01 E=2.3612e-09
これらの問題に対する万能の解決策というものはないのですが、例えば積分区間を分割することで改善することがあります。以下の例では積分区間を分割し、quadgk
では積分値の精度を、quadde
では評価回数を改善しています。
let f(x) = 1/cosh((10x-2))^2+1/cosh((100x-40))^4+1/cosh((1000x-600))^6
@printeval quadde(f, 0.0, 0.3, 0.5, 0.7, 1.0)
@printeval quadgk(f, 0.0, 0.3, 0.5, 0.7, 1.0)
end
;
quadde(f, 0.0, 0.3, 0.5, 0.7, 1.0): Evaluation=6976 I=2.108027355006400e-01 E=1.6566e-12 quadgk(f, 0.0, 0.3, 0.5, 0.7, 1.0): Evaluation=690 I=2.108027355005494e-01 E=1.9724e-09
次は $x=0$ 付近で急激に変化する関数です。
$$ \int_{0}^{\infty} \left( K_0(x) \right)^2 dx = \frac{\pi^2}{4} = 2.46740110027233\dots $$quadgk
は評価回数は多いものの、7桁目まで一致、9桁目を四捨五入すると8桁目まで正答しています。対して、quadde
は評価回数は少ないものの、 4桁まで一致と x = 0
付近の急な変化を見落としているようです。
let f(x) = besselk(0, x)^2
@printeval quadde(f, 0.0, Inf)
@printeval quadgk(f, 0.0, Inf)
plot(f, label=L"\left( K_0(x) \right)^2", xlim=[0.0, 2.0],
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, Inf): Evaluation=47 I=2.467369673162837e+00 E=1.9604e-08 quadgk(f, 0.0, Inf): Evaluation=765 I=2.467401095296906e+00 E=2.7633e-08
これも、やはり積分区間を区切ることで改善します。もとの積分区間 (0 ~ ∞) よりも狭い範囲 (0 ~ 2) を具体的に指定することで急な変化を見つけることができたようです。
let f(x) = besselk(0, x)^2
@printeval quadde(f, 0.0, 2.0, Inf)
@printeval quadgk(f, 0.0, 2.0, Inf)
end
;
quadde(f, 0.0, 2.0, Inf): Evaluation=79 I=2.467401098347747e+00 E=1.7313e-11 quadgk(f, 0.0, 2.0, Inf): Evaluation=750 I=2.467401093715945e+00 E=3.3688e-08
最後に積分区間に不連続な点を含む関数です。これは数値積分法にとっては悪夢のような問題といえます。
$$ \int_{0}^{1} \lfloor \min \left(\frac{x}{0.3}, 1\right) \rfloor dx $$とにかく、不連続に変化している $x = 0.3$ を発見できるかによって精度が大きく左右されます。この一点を見つけられなければまともな精度は得られません。
quadgk
は数百程度の評価回数で8桁を一致させています。quadde
は苦戦しており数万回の評価を経てなお4,5桁程度の一致にとどまります。
let f(x) = floor(min(x/3*10, one(x)))
@printeval quadde(f, 0.0, 1.0)
@printeval quadgk(f, 0.0, 1.0)
plot(f, label=L"\lfloor \min \left(\frac{x}{0.3}, 1\right) \rfloor",
xlim=[0.0, 1.0], xticks=0.0:0.1:1.0,
fillrange=0, fillalpha=0.1, fillcolor=:blue)
end
quadde(f, 0.0, 1.0): Evaluation=25787 I=6.999612521325285e-01 E=2.7833e-07 quadgk(f, 0.0, 1.0): Evaluation=675 I=7.000000081693309e-01 E=9.2649e-09
$x = 0.3$ で不連続になることが分かっていれば、ここを端点とすることで問題を回避できます。
let f(x) = floor(min(x/3*10, one(x)))
@printeval quadde(f, 0.0, 0.3, 1.0)
@printeval quadgk(f, 0.0, 0.3, 1.0)
end
;
quadde(f, 0.0, 0.3, 1.0): Evaluation=38 I=7.000000000000256e-01 E=3.1603e-10 quadgk(f, 0.0, 0.3, 1.0): Evaluation=30 I=7.000000000000000e-01 E=0.0000e+00