最終更新: 2015/11/8

特殊関数の精度保証

柏木 雅英

1. はじめに

本ライブラリが持つ精度保証付き数値積分solverや 精度保証付き常微分方程式solverを用いて、 特殊関数の精度保証を行う関数をいくつか試験的に作成している。

内部で数値積分や常微分方程式を解いているため、 基本的に遅い。 しかし、単純な原理に基づいているため精度保証されていることが 明快に理解でき、また計算に使う型を変えることによって高精度計算を 行うこともできるので、役に立つ場面もあるのではないかと考えている。 現状複素数に対応していない。 順次対応させるつもりではあるが、予定は未定である。

2. Gamma関数

gamma.hppの中に、 の実装がある。 いずれもinterval<T>を入力し、interval<T>を出力する。

詳細は、 ガンマ関数の精度保証付き計算メモ(Version: 2015/01/30) に書いた。

なお、gamma関数の極値点を計算するのにdigammaに対するKrawczyk法を行う関係で、 呼び出しの初回は非常に時間がかかる。

3. Lobachevsky関数

Λ(θ) = - ∫0θ log|2sin(t)|dt

で定義される関数をLobachevsky関数という。 lobachevsky.hppの中に、関数lobachevskyがある。 interval<T>を入力し、interval<T>を出力する。

詳細はまだ書いてない。

4. Airy関数

airy.hppの中に、 の実装がある。それぞれ、Ai関数、Bi関数、Ai関数の微分、Bi関数の微分である。 いずれもinterval<T>を入力し、interval<T>を出力する。

アルゴリズムは、単純にAiry関数を定義する初期値問題

f''(t) = t f(t)
を初期値 を使って解くだけ。 2階に直して解くので副産物として微分も得られる。 初期値問題を時刻0から単に解いていくだけなので、 絶対値の大きい入力に対しては精度も出ないし計算時間もかかる。

5. Bessel関数

bessel.hppの中に、第1種Bessel関数の実装がある。整数次のみ 非整数次も試験的に実装。

アルゴリズムは、積分 Jn(x) = 1/π * ∫0π cos(x sin t - nt) dt を計算するだけ。
非整数次の場合は、 Jnu(x) = (x/2)nu / sqrt(π) / gamma(nu + 0.5) * ∫0π cos(x * cos(t)) * sin(t)2nu dt を計算。nu>-0.5が必要。

6. Beta関数

beta.hppにある。引数を計算しやすいところにreduceした後、直接数値積分している。

7. ガウスの超幾何関数

hypergeom.hppにある。 ガウスの超幾何関数 2F1(a,b,c;z) を計算する。

オイラーの積分表示を直接計算している。そのため、 c > b > 0, |z| < 1 のようなパラメータの制限がある。

計算方法は hypergeom.pdf (Version: 2015/11/8) の通り。 途中に出てくる区間公比を持つ等比級数の和の公式は、 geoseries.pdf (Version: 2015/11/8) にある。