積分で定義される関数の高精度数値計算法：量子力学の数値解析から数学・数値解析へのフィードバック

Abstract

We applied a method for calculating indefinite integrals developed for atomic structure calculations to mathematical functions defined by integrals and obtained highly accurate results. The method is simple and has wide applicability. Examples of application are, elementary transcendental functions such as exponential, logarithmic, and arctangent functions, then error function, Fresnel integrals, exponential, cosine, and sine integrals, and incomplete gamma functions.

1 序論

我々は量子力学の一分野である原子構造計算につき，高精度計算を実現するため，既存の数値計算法を利用するだけでなく，新しい計算法を考案し，課題を克服してきた [1]．ここでは，原子構造計算で開発した不定積分の計算法を一般の数学関数の計算に適用することにより，積分で表わされる関数を高精度で数値計算する一般的な方法を報告する．

2 数学的基礎

積分で表わされる関数の数学的基礎は，解析学でよく知られているように，微分積分学の基本定理である　${\int }_a^{b}f\left(x\right)dx=F\left(b\right)-F\left(a\right)$．ここで左辺は区間$\left[a,b\right]$における関数$f\left(x\right)$のRiemann積分である．また，右辺は原始関数$F\left(x\right)$の両端での値の差であり，原始関数$F\left(x\right)$と導関数$f\left(x\right)$との関係は$(d/dx)F\left(x\right)=f\left(x\right)$である [2,3,4]．

この定理で，積分区間の上限$b$を変数$x$とすると，$F\left(x\right)=F\left(a\right)+{\int }_a^{x}f\left(x\right)dx.$　或いは，積分区間の下限$a$を変数$x$とすると，$F\left(x\right)=F\left(b\right)-{\int }_x^{b}f\left(x\right)dx.$　これらの$F\left(x\right)$を数値計算する一般的な方法を以下に述べる．これは$F\left(x\right)$が初等関数で表わせない場合にも対応できる．更に，$f\left(x\right)$が数値データでのみ与えられる場合にも拡張して使える．

3 計算方法

関数が定義されている全区間を等間隔で$N$分割する，$x_j,$$j=0,\mathrm{...}N.$その一つの小さな区間$(x_k,x_{k+1})$ における積分 ${\int }_{x_k}^{x_{k+1}}f\left(x\right)dx$ (単一区間の積分と呼ぶことにする)を考える．この積分をメッシュ点での関数値のみを使って計算する [1]．隣接する$(n+1)$個の分点をとる．これらの分点での関数値を用いてLagrange補間多項式を構成する．単一区間の積分公式はLagrange補間多項式を単一区間で積分することにより得られる．得られた積分公式では積分区間の外側の分点も使う．(積分公式の次数は対応する補間多項式の次数と同じ$n$である．)特に，中央の区間でのLagrange補間は高精度であるから [5,6]，その一区間での積分も高精度である [1]．この積分は非常に有用なので，それを特に中央区間の積分公式(central-interval integration formula, CIIF)と呼ぶことにする [1]．更に，単一区間の積分値の累積を計算することにより，メッシュ点$x_k$における積分値を計算する   

$F\left(x_k\right)=F\left(a\right)+{\int }_a^{x_k}f\left(x\right)dx=F\left(a\right)+\sum _{j=1}^k{\int }_{x_{j-1}}^{x_j}f\left(x\right)dx,$

或いは   

$F\left(x_k\right)=F\left(b\right)-{\int }_{x_k}^{b}f\left(x\right)dx=F\left(b\right)-\sum _{j=k}^{N-1}{\int }_{x_j}^{x_{j+1}}f\left(x\right)dx.$

ここで$b$は$N$番目のメッシュ点とする．メッシュ点の中間の点における$F\left(x\right)$は補間により計算する．(高精度の多項式補間法は [5,6]に示されている．)

更に，関数が区間$\left[0,\infty \right)$で定義されている場合，原点近傍では冪級数展開式を，無限遠点近傍では漸近級数展開式により$F\left(x\right)$を計算する．それらの中間の領域ではCIIFによる積分とその累積和を計算する．計算の際のパラメータは，分割の区間幅，数値積分公式の次数，冪級数展開を使う領域の上限の位置，漸近展開を使う領域の下限の位置，である．これらは要求精度に応じてテストをしながら決める．

4 計算例

計算した関数は以下の通りである．関数の定義は [7,8]に従っている．

(a) 積分が初等関数で表わせる，或いは初等関数となるもの

指数関数，対数関数，逆三角関数(特に逆正接関数)これらは，コンパイラの組み込み関数を利用して，計算精度の評価を行った．

(b) 初等関数で表わせない積分

誤差関数と補誤差関数，Fresnel積分，積分指数関数，積分余弦関数，積分正弦関数，不完全ガンマ関数

これらの関数は他の計算法によるプログラム [9]を利用して，計算精度の評価を行った．

5 計算結果

典型例として対数関数を示す．対数函数は次式で定義される   

$F\left(x\right)=\ln \left(x\right)={\int }_1^x(1/t)dt.$

計算は倍精度演算で行った．$x$を変えながら，CIIF積分とその累積和を計算した．CIIF積分の相対誤差(の絶対値)をFigure 1に示す．区間幅$h$をパラメータとして，公式の次数を上げると誤差は減少し，倍精度目一杯の精度が出る．(尚，相対誤差が小さな所で，誤差が下がらなくなるのは，丸め誤差の影響による．)同じ次数で$h$を小さくすると，誤差は減少する．

Figure 1.

 Relative errors of integrals for single interval by using central-interval integration formula (CIIF) vs. degree of the integration formula. $h$ is the width of the interval. Interval for the integration is $[2-h,2]$.

積分の累積和，ここでは$x=2$での関数値$\ln \left(2\right)$，の相対誤差をFigure 2に示す．積分公式の次数を上げると誤差は減少し，倍精度目一杯の精度が出る．($h=1/64$で積分公式の次数9では，積分の累積和の計算値が組み込み関数の計算結果と全桁で一致し，相対誤差はゼロになった．縦軸を対数スケールで表示すると，ゼロの対数はマイナス無限大となって不都合が生じるが，使用ソフト(KaleidaGraph Version 4.5)の仕様では縦軸の範囲を図のように指定すると，下に急降下する直線で表示される．)

Figure 2.

 Relative errors of $\ln \left(2\right)$ that is a cumulative sum of integrals for single interval by using central-interval integration formula (CIIF) vs. degree of the integration formula.

他の関数についても，同様の計算精度を得ている．関数の性質，例えば特異性や急激な変化，に応じた工夫が必要な場合には，対処している．これらについては別途述べる．

6 議論

既存の方法との関係を述べる．関数計算は過去数世紀にわたって調べられてきた．総合報告には，例えば， [7,10]がある．今では種々の計算法がある．典型的には全区間をいくつかの区間に分けて，それぞれの区間で適切な計算法を用いる．計算法の分類と特徴を以下に示す：

◇原点近傍

・冪級数展開を使う

◇無限遠点近傍

・漸近展開を使う

◇中間の領域

・連分数展開を使う

・漸近展開の領域から更に内側へ計算可能な領域を広げる

・数値積分を使う

・元の被積分関数をそのまま使う

・Romberg積分，Gauss型公式，等を使う

・単一区間の積分とその累積を計算する

・積分の表式を変換してから積分する

・無限区間の積分に変形して台形公式を使う

・多項式近似(全区間をいくつかの区間に分割し，その区間で近似式を作る)を使う

・最大誤差を最小にする近似法(Chebyshev近似)を使う

・Taylor級数展開法(展開係数(高階微分)を数値で求めて表で持つ)を使う

・有理多項式近似を使う

冪級数展開と漸近展開は確立された方法であり，我々も利用している．新しい方法が出てくるのは中間の領域である．今回の計算法は，単一区間の積分とその累積を計算する，という新しいものである．単一区間の積分法は単純かつ高精度という特徴を持つ．これは適用性が広いというメリットがある．

参考文献

• 1   H. Ishikawa, Numerical methods of calculation for atomic structure calculations, unpublished.
• 2   T. Takagi, Kaiseki Gairon (Introduction to Analysis), Iwanami, Tokyo, 1961.
• 3   V. I. Smirnov, Koto Sugaku Kyotei, Vol. 1, Kyoritsu, Tokyo, 1958. (A Course of Higher Mathematics, Vol. 1, Pergamon, Oxford, 1964.)
• 4   E. Hairer, G. Wanner, Kaiseki Kyotei, Springer Japan, Tokyo, 1997. (Analysis by Its History, Springer, New York, 1996.)
• 5   H. Ishikawa, J. Phys. A, 35, 4453 (2002).
• 6   H. Ishikawa, J. Comput. Chem. Jpn., 6, 199 (2007).
• 7   M. Abramowitz, I. A. Stegun, ed., Handbook of Mathematical Functions, Dover, New York, 1972.
• 8   S. Moriguti, K. Udagawa, S. Hitotumatu, Sugaku Koshiki III (Mathematical Formulas III), Iwanami, Tokyo, 1960.
• 9   T. Watanabe, M. Natori, T. Oguni, ed., Fortran77 Niyoru Suti Keisan Sofutowea (Numerical Calculation Software with Fortran77), Maruzen, Tokyo, 1989.
• 10   F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark, NIST Handbook of Mathematical Functions, Cambridge University Press, 2010.