数値計算

2008.11.4 平野拓一

1. 数値計算とは

コンピュータを利用して数値的に値を評価する手法。例えば、不定積分が求まらない場合に区分求積で定積分値を求めるなどである。数値計算はコンピュータを利用するが、一般に、

精度は高い方が良い
計算時間は高速な方が良い

のは誰でも望むことである。しかし、上記の２つは相反する要求となる。上記２つの要求を満たすようなアルゴリズムを開発するのが数値計算の専門家の仕事となる。

2. 補間と補外

数値的に得られるデータは離散的である。その間の値、あるいはその外の値を知るためにサンプリングした離散値から連続関数を生成することがしばしば行われる。そして、生成した連続関数から間の値を予測することを補間(interpolation)、サンプリング範囲外の値を予測することを補外(extrapolation)と言う。

2.1 ラグランジュの多項式、ラグランジュ補間（多項式近似）

n点ある場合はn-1次の多項式で近似する方法。

2.1.1 スプライン補間（3次のラグランジュ補間）

(作成中）

2.2 Pade近似（有理関数近似）

有理関数で近似する方法。ラグランジュ近似では表現不可能な、極（n次の極とは、無限大に発散するが1/x^nで表すことができる特異点）がある関数を近似するのに精度が良い。

上式の未知数の数はp0,…pμのμ+1個と、q1,…qνのν個(q0は定数項のp0/q0を表すためなので、任意[例えば1]に選べる)で、合計μ+ν+1個である。よって、サンプル点がn個の場合はn=μ+ν+1と選べばよい。上式の係数は次のように求めることができる。

3. 数値積分

まず、1次元積分について説明する。

3.1 等間隔分割の場合

3.1.1 短冊和（リーマン積分の定義）

長方形の面積の和。

[参考文献]
・num_diff_int.pdf （数値微分と数値積分 by T. Hirano）

3.1.2 台計公式（直線補間を利用）

台形面積の和。

[参考文献]
・num_diff_int.pdf （数値微分と数値積分 by T. Hirano）

3.1.3 シンプソンの公式（２次曲線補間を利用）

３点を用い、２次曲線で補間して面積を計算する。滑らかな関数で精度が高い。

[参考文献]
・num_diff_int.pdf （数値微分と数値積分 by T. Hirano）

3.2 不等間隔分割の場合

3.2.1 ガウス・ルジャンドル積分（１次元の場合）

ルジャンドルの多項式の直交性を利用した不等間隔数値積分。被積分関数が2n-1次の多項式で表現されるならば、分点数をn点以上にすると誤差は0となる。有限要素法では滑らかな関数を積分するので、よく用いられる。

[アルゴリズム]

を求めたいとする。積分範囲は-1～1であるが、最後に積分範囲の変換式を示すように、これでも一般性を失わない。f(x)を2n-1次の多項式とすると、n次のルジャンドルの多項式P_n(x)で割った商Q(x)と剰余R(x)で表す。

ところで、ルジャンドルの多項式は互いに直交している。

なぜならば、Q(x)はn-1次であるから、であり、直交性から上式の下線の積分は0となる。

求積を行う際、下記のように等式で結べるとする（f(x)は2n-1次の多項式であると仮定しているから）。

Iを求めるためにはR(x)を求めなければならないように見えるが、より、Pn(x)の０点をサンプリング点に選ぶととなる（R(x)を求める必要がない）。

さて、残るはwiはどのような値であるか定めることである。上の右２式より、

R(x)は任意のn-1次の多項式に対して成り立たなければならない。R(x)として、P₀, P₁, ..., P_n-1のn個の関数を選ぶと次のn個の連立一次方程式が得られる。

w₁, ..., w_nは上式を解くと求めることができる。これでn次の公式が得られた。

例として、n=3の場合のx_iとw_iを示す。

積分範囲をa～bにするには次のように座標変換するだけでよい。

ただし、

方形領域であれば、多次元への拡張も容易である。

[参考文献]

num_int_gauss_legendre.pdf (Mathematicaによる確認 by T. Hirano)
G.R. Cowper, "Gaussian Quadrature Formulas for Triangles," Int. J. Numer. Mathods Eng., Vol.7, pp.405-408, 1973. （三角形領域に対称にサンプリング点を配置した場合のガウス・ルジャンドル求積法）
以下、私（平野）の理解で論文を要約する。三角領域Sでの面積分を次のように表すことができるとする。
(※)
ここで、は三角領域の面積である。論文の目的は対称な位置にサンプリング点を配置することなので、面積座標で書いた場合、の形のものは

の６点に配置しなければならない。同様にの形のものは

の３点に配置しなければならない。全ての面積座標が同じ場合は中心の１点である。このように、座標は３点あるいは６点ペアで追加されることとなる。さて、問題は式(※)の右辺の係数w_i,ξ_i,η_i(面積座標の条件より、ξ_i,η_i,ζ_iのうち1文字は従属である)およびサンプリング点などをどのように決定するかである。
例えば、およびの７点をサンプリング点を用いる公式を得る場合には次のようにする。未知数の数はw_i,ξ_i,η_iおよびα,β,γの6文字である。の6個の関数は全て式(※)を満たさなければならないので、これらを式(※)に代入して6元連立方程式を得る。もちろん、これは非線形で解きにくいが、本論文によると何とかこれを解くことができ、サンプリング点と重みを決定することができる（そうである）。
H.T. Rathod, K.V. Nagaraja, and B. Venkatesudu, "Symmetric Gauss Legendre quadrature formulas for composite numerical integration over a triangular surface," Elsevier, Applied Mathematics and Computation, Vol.188, No.1, pp.865-876, 1 May 2007. （三角形領域を正方形領域に座標変換してガウス・ルジャンドル求積法を適用する方法）
P.C. Hammer, O.P. Marlowe, and A.H. Stroud, "Numerical Integration Over Simplexes and Cones," Math. Tables Aids Comp. 10, 130-7, 1956. （四面体のサンプリング点と重み）