Julia 言語で汎用一次元数値積分関数を書きました

二重指数関数型数値積分公式 を使った、非適応型の数値積分プログラムです。

GitHub - machakann/DoubleExponentialFormulas.jl: One-dimensional numerical integration using the double exponential formula

高橋・森によって提案された変数変換を施すことによって、被積分関数を積分区間の両端で急激に減衰する新しい関数に変換します。このような関数は基本的に元の関数よりも数値積分が容易で、主に台形則などを使って高速かつ精度よく近似解を得ることができます。

被積分関数が積分区間全体にわたってなめらかに変化する場合に、特に少ない評価点数で高精度な値を返します。また、もう一つの特徴として端点での特異性に非常に強い点が挙げられます。

Tanh-sinh quadrature

積分区間や被積分関数の振る舞いによって多くの変形がありますが、最も汎用的かつ基本的な三種を実装しました。tanh-sinh quadrature はその一つで、任意の被積分関数の [-1, 1] の区間での数値積分を求めます。

\begin{align} \int_{-1}^{1} f(x) dx = \int_{-\infty}^{\infty} f\left(x(t)\right) \frac{dx}{dt} dt = \int_{-\infty}^{\infty} f\left(x(t)\right)w(t) dt \end{align}

 x についての被積分関数  f(x) を、  t についての新しい被積分関数  f(x(t))w(t) へ変換し、それに伴い積分区間をうつしています。ここで、 x(t) および  w(t) は以下のような関数です。

\begin{eqnarray} x\left(t\right) = \sinh \left( \frac{π}{2} \sinh t \right) \\ w\left(t\right) = \frac{π}{2} \cosh t \cosh\left(\frac{π}{2} \sinh t \right) \end{eqnarray}

この変数変換が肝なので、以下の積分を例にとり被積分関数が実際にどのように変化するかを見てみましょう。

\begin{align} \int_{-1}^{1} \frac{dx}{\left(2-x\right)\left(1-x\right)^{1/4}\left(1+x\right)^{3/4}} \end{align}

この関数は積分区間の両端  x = \pm1 で急速に発散しており、数値積分を難しくしています。しかし、変数変換を施したあとの関数  f(x(t))w(t) は正負の両方向で急速に減衰する関数に変換されていることがわかります。

fig.1

元の関数の区間両端での発散が解消された代わりに両無限区間 [-∞, ∞] での積分になっています。しかし、変換後の被積分関数は正負の両方向へ向かうにつれ急激に減衰する関数となっているため、台形則の適用により無限和に近似し  f(x_k)w_k が十分小さくなった地点で打ち切ることができます。

\begin{align} \int_{-\infty}^{\infty} f\left(x(t)\right)w(t) dt \approx h\sum_{k = -\infty}^{\infty} f(x_k)w_k \end{align}

ここで  x_k および  w_k は台形則における刻み幅  h を用いて次のように表せます。

\begin{eqnarray} x_k = x(kh) \\ w_k = w(kh) \end{eqnarray}

 x_k および  w_k は被積分関数に依存しないため、刻み幅  h を決めておけば事前に  x_k,  w_k のテーブルを計算しておくことができます。これも速度が求められる場面では重要な特性といえるでしょう。また事前に計算しない場合でも、tanh-sinh quadrature については  x(t) が奇関数、 w(t) が偶関数なので、 t = 0 に対する対称性 ( x(t) = -x(-t), w(t) = w(-t)) を利用して計算量を減らすことが可能です。

さらに別の変数変換を適用することで任意の有限区間 [a, b] での積分をすべてカバーします。

\begin{align} \int_{a}^{b} f(x) dx = \frac{b - a}{2} \int_{-1}^{1} f\left(x(u)\right) du = \frac{b - a}{2} \int_{-\infty}^{\infty} f\left(x\left(u(t)\right)\right)w(t) dt \end{align}

\begin{align} where\ \ x(u) = \frac{b + a}{2} + \frac{b - a}{2}u \end{align}

DoubleExponentialFormulas.jl はほかに exp-sinh quadrature, sinh-sinh quadrature を組み合わせて任意の有限区間 [a, b]、片無限区間 [a, ∞] [-∞, b]、両無限区間 [-∞, ∞] に対応しています。

使い方

関数 quaddeFloat64 の精度で数値積分を計算します。

    I, E = quadde(f::Function, a::Real, b::Real, c::Real...;
                  atol::Real=zero(Float64),
                  rtol::Real=atol>0 ? zero(Float64) : sqrt(eps(Float64)))

I は得られた積分値、E は推定誤差です。たとえば、f(x) = 1/(1 + x²) の区間 [-1, 1] での積分は次のように計算します。

using DoubleExponentialFormulas
using LinearAlgebra: norm

f(x) = 1/(1 + x^2)
I, E = quadde(f, -1.0, 1.0)

I ≈ π/2                   # true
E ≤ sqrt(eps(I))*norm(I)  # true

atolrtol はそれぞれ目標とする絶対誤差および相対誤差で、E がこれらのうち大きいほうを下回るまで精度を上げながら反復を繰り返します。ただし、E はあくまで推定値であり I と真の値の正確な誤差ではありません。しかし、少なくとも I が収束した場合には E <= max(atol, rtol*norm(I)) が真となることが期待できます。言い換えれば、その条件が満たされていない場合、I は既定の反復回数のうちに収束しなかったことを意味しており、信頼できない値であると考えられます。

基本的には積分値の大きさに応じた精度を要求できる rtol を使うのが便利です。しかし、予想される積分値 I の絶対値が非常に小さい場合には rtol*norm(I) がアンダーフローを起こす恐れがあります。こうなると収束条件を満たせず反復回数が大きくなってしまうので、これを避けるために同時に atol にも小さい値を設定しておくと計算量を節約できる場合があります。atol が与えられた場合は、rtol は明示的に与えられない限り使われないので両方を明示的に指定しましょう。

区間には Inf を指定することができます。

# Computes ∫ 1/(1+x^2) dx in [0, ∞)
I, E = quadde(x -> 1/(1 + x^2), 0.0, Inf)
I ≈ π/2    # true

# Computes ∫ 1/(1+x^2) dx in (-∞, ∞)
I, E = quadde(x -> 1/(1 + x^2), -Inf, Inf)
I ≈ π      # true

fig.2

また、a, b 以降にも区間が与えられた場合には積分区間を分割して計算した和を返します。

\begin{align} I = \int_{a}^{b} f(x) dx + \int_{b}^{c} f(x) dx + \cdots \end{align}

これは、積分区間に不連続点・発散する点を含む場合に有用です。例えば、関数  {f(x) = \frac{1}{\sqrt{|x|}}} {x = 0} で発散しますが、区間を分割し  {x = 0} を端点にすることでうまく計算できるようになります。

# Computes ∫ 1/sqrt(|x|) dx in (-1, 1)
# The integrand has a singular point at x = 0
I, E = quadde(x -> 1/sqrt(abs(x)), -1.0, 0.0, 1.0)
I ≈ 4    # true

fig.3

得意・不得意な関数の傾向

基本的には被積分関数  {f(x)} が積分区間全域にわたって緩やかに変化するような関数は効率よく計算するようです。逆に積分区間の一部だけ激しく変化するような関数の場合は、最も変化の激しい部分に引っ張られて全体の評価点数を増やすので効率が良くないです。これは非適応型の宿命といえます。具体的な例を挙げると、積分区間にスパイクが発生するような関数は評価回数が多くなってしまいます。

fig.4

また、上でも述べられたように端点の特異性に関してはほかのアルゴリズムと比べても強いため、ほかのアルゴリズムでは対応できない場合に試してみるのはよいかもしれません。

QuadGK との比較

→nbviewerでみる

gist.github.com

感想など

書く前は難しそうだと思っていたので、思ったよりもちゃんと動くな、というのが正直な感想です。いままで数値積分関数はブラックボックスとして扱っていたので勉強になりました。数値積分関数の説明には往々にして注意書きがいろいろあったり、「必ずしも正しい値を返すとは限りません」みたいな歯切れの悪い文言が並んでいたりするものなんですけれども、なんとなく理由がわかったような気がします。上のスパイクが発生するような関数なんかは二重指数関数型数値積分法に限らず数値積分関数にとっては悪意の塊みたいな関数なんだな、と。

汎用性のために一部速度や精度を犠牲にしている部分もあります。また、二重指数関数型積分法には被積分関数によっていろいろな変種が考案されているので、汎用の関数を書いておいてなんですが問題にあわせて最適化したものをつかうのが一番ですね。それでもほとんどの被積分関数に関して7~8桁(デフォルト)の精度での計算には問題なく動くようです。13桁かそれ以上を目指す場合には被積分関数によっては厳しいかもしれないです。

今のところはかなり素朴な実装 (v0.0.5) になっていて特にひねりもないので、何か思いついたら改良していこうと思います。