Gauss-Legendre 積分公式

 


Gauss-Legendre 積分公式

Gauss-Legendreの積分公式は,有限要素法のアイソパラメトリック要素における数値積分などに用いられており,積分点および重みを定義すれば計算フローが単純で精度も比較的高い.

Gaussの積分公式では,積分範囲は $-1 \leqq t \leqq 1$ の範囲で定義され,適切な積分点 $t_i$ と重み $w_i$ を選定することにより,以下のように表現される.

\begin{equation} \int_{-1}^{1} g(t) dt = \sum_{i=1}^n w_i g(t_i) \end{equation}

上式を任意の積分範囲 $a \leqq x \leqq b$ に適用するため,以下に示す積分の変数変換を用いる.

\begin{equation} \int_a^b f(x) dx = \int_{-1}^{1} f(g(t))\cdot g'(t) dt = \cfrac{b-a}{2} \int_{-1}^{1} f\left(\cfrac{b-a}{2} t + \cfrac{a+b}{2}\right) dt \qquad\qquad g(t)=\cfrac{b-a}{2} t + \cfrac{a+b}{2} \end{equation}

上の変数変換された表現に,Gaussの積分公式を適用するとともに,Legendreの多項式により積分点と重みを決定することとすれば,以下に示す Gauss-Legendreの積分公式が得られる.

\begin{align} & \int_a^b f(x) dx \doteqdot \frac{b-a}{2}\sum_{i=1}^n w_i\cdot f(x_i) \\ & x_i=\frac{b-a}{2}\cdot t_i + \frac{a+b}{2} \\ & w_i=\frac{2(1-t_i{}^2)}{n^2 P_{n-1}{}^2 (t_i)} \\ & t_i : \text{$P_n(t_i)=0$ を満足する実数} \end{align}

ここで, $P_{n}(x)$ はLegendreの多項式であり,以下のとおり定義される.

\begin{equation} P_{n}(x)=\frac{1}{2^n n!}\frac{d^n}{dx^n}(x^2-1)^n \qquad (-1\leqq x \leqq 1) \end{equation}

上式およびその微分は以下の漸化式にて表現できる.

\begin{align} &P_{n}(x)=\frac{(2 n-1) x P_{n-1}(x)-(n-1) P_{n-2}(x)}{n} \qquad P_0(x) =1 \qquad P_1(x)=x \\ &P'_{n}(x)=\frac{n \left\{x P_{n}(x)-P_{n-1}(x)\right\}}{x^2-1} \end{align}

実際にこの積分を実行するには,まず積分点数 $n$ ,積分点 $t_i$ ( $P_n(t_i)=0$ を満足する $t_i$ ) および重み $w_i$ を求めなければならない.

$P_n(t_i)=0$ を満足する $t_i$ は, $P_n(x)$ およびその導関数を漸化式より求められるため,Newton法により求めることができる.

\begin{equation} t_{i(k+1)}=t_{i(k)}-\frac{P_n(t_{i(k)})}{P'_n(t_{i(k)})} \end{equation}

ここで, $t_i$ の値は,正数側と負数側で対称となるため,正数側の値のみを求め,符号を変更して $n$ 個の積分点の値を求めればよい. この際, $n$ が奇数の場合は $t_i=0$ を加えることを忘れないようにする.

Newton法適用時の, $t_i$ の初期値としては, $t_{i(0)}=\cos\left(\cfrac{\pi}{2}\cfrac{4 i-1}{2 n+1}\right)$ を用いればよい.

以下に各 $n$ に対応する積分点 $t_i$ と重み $w_i$ の計算事例を示す.

n i ti wi
2 1 -0.577350269189626D0 1.000000000000000D0
2  0.577350269189626D0 1.000000000000000D0
3 1 -0.774596669241483D0 0.555555555555555D0
2  0.000000000000000D0 0.888888888888889D0
3  0.774596669241483D0 0.555555555555555D0
4 1 -0.861136311594053D0 0.347854845137454D0
2 -0.339981043584856D0 0.652145154862546D0
3  0.339981043584856D0 0.652145154862546D0
4  0.861136311594053D0 0.347854845137454D0
5 1 -0.906179845938664D0 0.236926885056189D0
2 -0.538469310105683D0 0.478628670499366D0
3  0.000000000000000D0 0.568888888888889D0
4  0.538469310105683D0 0.478628670499366D0
5  0.906179845938664D0 0.236926885056189D0
10 1 -0.973906528517172D0 0.066671344308688D0
2 -0.865063366688985D0 0.149451349150580D0
3 -0.679409568299024D0 0.219086362515982D0
4 -0.433395394129247D0 0.269266719309996D0
5 -0.148874338981631D0 0.295524224714753D0
6  0.148874338981631D0 0.295524224714753D0
7  0.433395394129247D0 0.269266719309996D0
8  0.679409568299024D0 0.219086362515982D0
9  0.865063366688985D0 0.149451349150580D0
10  0.973906528517172D0 0.066671344308688D0


Gauss-Legendre 積分公式を用いた完全楕円積分

Gauss-Legendre 積分公式を用いて第一種および第二種完全楕円積分の数値解を計算した事例を示す.

ここに,第一種完全楕円積分 $K(p)$ および第二種完全楕円積分 $E(p)$ は以下のとおり定義される.

\begin{equation} K(p)=\int_0^{\pi/2}\cfrac{d\theta}{\sqrt{1-p^2 \sin^2\theta}} \qquad E(p)=\int_0^{\pi/2}\sqrt{1+p^2 \sin^2\theta}d\theta \end{equation}

以下の事例に示したとおり,積分点数が少なくても,積分区間を細かく分割することにより,精度を高めることができる.

Gauss-Legendre 積分による結果 (積分点 = 2)

積分点数を 2 として,積分区間 $0\sim 2\pi$ を細かく分割していき,積分値の変化が $10^{-10}$ 以下となったところで計算を終了した事例. 条件の厳しい $K(0.9999)$ のケースにおいては,323 分割で,指定した精度となり計算を終了している.

Gauss-Legendre rule (n=2)
ita. p K(p) E(p)
2 0.0000000E+000.1570796E+010.1570796E+01
3 0.1000000E+000.1574746E+010.1566862E+01
4 0.2000000E+000.1586868E+010.1554969E+01
4 0.3000000E+000.1608049E+010.1534833E+01
5 0.4000000E+000.1640000E+010.1505942E+01
6 0.5000000E+000.1685750E+010.1467462E+01
6 0.6000000E+000.1750754E+010.1418083E+01
7 0.7000000E+000.1845694E+010.1355661E+01
9 0.8000000E+000.1995303E+010.1276350E+01
13 0.9000000E+000.2280549E+010.1171697E+01
323 0.9999000E+000.5645148E+010.1000515E+01
time= 0.12500000 (sec)

Gauss-Legendre 積分による結果 (積分点 = 10)

積分点数を 10 として,積分区間 $0\sim 2\pi$ を細かく分割していき,積分値の変化が $10^{-10}$ 以下となったところで計算を終了した事例. 条件の厳しい $K(0.9999)$ のケースでも,53 分割で,指定した精度となり計算を終了している.

Gauss-Legendre rule (n=10)
ita. p K(p) E(p)
2 0.0000000E+000.1570796E+010.1570796E+01
2 0.1000000E+000.1574746E+010.1566862E+01
2 0.2000000E+000.1586868E+010.1554969E+01
2 0.3000000E+000.1608049E+010.1534833E+01
2 0.4000000E+000.1640000E+010.1505942E+01
2 0.5000000E+000.1685750E+010.1467462E+01
2 0.6000000E+000.1750754E+010.1418083E+01
2 0.7000000E+000.1845694E+010.1355661E+01
3 0.8000000E+000.1995303E+010.1276350E+01
3 0.9000000E+000.2280549E+010.1171697E+01
53 0.9999000E+000.5645148E+010.1000515E+01
time= 0.0000000 (sec)


Simpson 公式との比較

Gauss-Legendre 積分と Simpson 公式による積分値の比較を行った.計算内容は以下の式により円周率を計算すること.

\begin{equation} \pi=2\int_{-1}^{1}\sqrt{1-x^2}dx \end{equation}

Simpson 公式は以下の型式を使用した.以下は分割された比較的小さな1要素に対するものとすると,要素の始点・終点および中間点の3点の関数値により積分値を評価しているとみなすことができる. 実計算では,逐次下式で計算された値を加算していき,最終的な積分値を求めている.

\begin{equation} \int_{a}^{b}f(x) dx \doteqdot \frac{b-a}{6}\left[f(a)+4f\left(\cfrac{a+b}{2}\right)+f(b)\right] \end{equation}

Gauss-Legendre 積分では,1要素あたりの積分点数は Simpson 公式との整合性(3点の値で小さな領域の積分値を評価する)を考慮して 3 とした.

計算結果は以下のとおり.n は積分範囲の分割数,Estimated value は積分値,Error は正解 ( 4*atan(1) の倍精度値と設定) に対する差分.

Simpson
n Estimated value Error
1 0.266666666667E+01 -0.474925986923E+00
10 0.312700815870E+01 -0.145844948866E-01
100 0.314113320534E+01 -0.459448250567E-03
1000 0.314157813021E+01 -0.145233758078E-04
10000 0.314159219438E+01 -0.459214426485E-06
100000 0.314159263907E+01 -0.145226861648E-07
1000000 0.314159265314E+01 -0.454499993197E-09

Gauss - Legendre
n Estimated value Error
1 0.318323451563E+01 0.416418620406E-01
10 0.314286934102E+01 0.127668742748E-02
100 0.314163289136E+01 0.402377679936E-04
1000 0.314159392559E+01 0.127200198508E-05
10000 0.314159269381E+01 0.402228939045E-07
100000 0.314159265486E+01 0.127201404965E-08
1000000 0.314159265363E+01 0.402087252382E-10


Programs

Programs

FilenameDescription
f90_GAUSSIP.txtprogram for calculation of Gauss points and weights for integration
f90_CEIG.txtprogram for calculation of Complete elliptic integration
f90_SIMPSON.txtprogram for calculation of PI by Simpson rule
f90_GINT3.txtprogram for calculation of PI by Gauss-Legendre rule with 3 Gauss points


pic
inserted by FC2 system