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=2, 3, 4, 5, 10$ に対応する積分点 $t_i$ と重み $w_i$ の計算プログラムと計算結果事例を示す．


FilenameDescription
py_gauss.pyGauss-Legendre Abscissas and Weights

上記プログラムは，理論に忠実に計算するものであるが，numpyには積分点と重みを計算する関数が準備されており，以下のコードでそれらを計算・表示できる．簡単！


# Gauss-Legendre Abscissas and Weights
import numpy as np
print('{0:>5s} {1:>15s} {2:>15s}'.format('i','t','w'))
for n in [2,3,4,5,10]:
    t,w=np.polynomial.legendre.leggauss(n)
    for i in range(0,n):
        print('{0:5d} {1:15.12f} {2:15.12f}'.format(i+1,t[i],w[i]))
    print()


   i               t               w
    1 -0.577350269190  1.000000000000
    2  0.577350269190  1.000000000000

    1 -0.774596669241  0.555555555556
    2  0.000000000000  0.888888888889
    3  0.774596669241  0.555555555556

    1 -0.861136311594  0.347854845137
    2 -0.339981043585  0.652145154863
    3  0.339981043585  0.652145154863
    4  0.861136311594  0.347854845137

    1 -0.906179845939  0.236926885056
    2 -0.538469310106  0.478628670499
    3  0.000000000000  0.568888888889
    4  0.538469310106  0.478628670499
    5  0.906179845939  0.236926885056

    1 -0.973906528517  0.066671344309
    2 -0.865063366689  0.149451349151
    3 -0.679409568299  0.219086362516
    4 -0.433395394129  0.269266719310
    5 -0.148874338982  0.295524224715
    6  0.148874338982  0.295524224715
    7  0.433395394129  0.269266719310
    8  0.679409568299  0.219086362516
    9  0.865063366689  0.149451349151
   10  0.973906528517  0.066671344309


２重積分

Gauss-Legendre積分の積分点および積分値は，積分領域が $-1 \leqq u \leqq 1$，$-1 \leqq v \leqq 1$ で与えられるため，任意の領域を積分するためには，変数変換を考える必要がある． 一般に，Jacobi行列を$[J]$として

\begin{equation} \int\int f(x,y)dx dy=\int\int f\left(x(u,v),y(u,v)\right) \det[J] du dv \end{equation}

ここで， $x_a \leqq x \leqq x_b$，$y_a \leqq y \leqq y_b$ の長方形領域を考えた場合，Gauss-legendre積分では

\begin{equation} x(u,v)=\cfrac{x_b-x_a}{2}u+\cfrac{x_a+x_b}{2} \qquad\qquad y(u,v)=\cfrac{y_b-y_a}{2}v+\cfrac{y_a+y_b}{2} \end{equation}
であるから，Jacobi行列 $[J]$ およびJacobian(jacobi行列の行列式) $\det[J]$ は，
\begin{equation} [J]= \begin{bmatrix} \cfrac{\partial x}{\partial u} & \cfrac{\partial x}{\partial v} \\ \cfrac{\partial y}{\partial u} & \cfrac{\partial y}{\partial v} \end{bmatrix} = \begin{bmatrix} \cfrac{x_b-x_a}{2} & 0\\ 0 & \cfrac{y_b-y_a}{2} \end{bmatrix} \qquad\qquad \det[J]=\cfrac{(x_b-x_a)(y_b-y_a)}{4} \end{equation}

となる． したがって，Gauss-Legendre積分の定義により

\begin{equation} \int_{x_a}^{x_b}\int_{y_a}^{y_b} f(x,y) dx dy \doteqdot \sum_{i=1}^n \sum_{j=1}^n \left\{w_i w_j f\left(x(u_i),y(v_j)\right)\right\} \cdot \cfrac{(x_b-x_a)(y_b-y_a)}{4} \end{equation}

$S'=\sum_{i=1}^n \sum_{j=1}^n \left\{w_i w_j f(x(u_i),y(v_j))\right\}$ と置き，４積分点モデル，９積分点モデルでの $S'$ の求め方を以下に示す．

４積分点
\begin{align} S'=&w_1 \left\{w_1 f(x(u_1),y(v_1)) + w_2 f(x(u_1),y(v_2))\right\}\\ +&w_2 \left\{w_1 f(x(u_2),y(v_1)) + w_2 f(x(u_2),y(v_2))\right\}\\ \end{align} 

i  j              wi              wj               ui              vj
1  1  1.000000000000  1.000000000000  -0.577350269190 -0.577350269190
1  2  1.000000000000  1.000000000000  -0.577350269190  0.577350269190
2  1  1.000000000000  1.000000000000   0.577350269190 -0.577350269190
2  2  1.000000000000  1.000000000000   0.577350269190  0.577350269190
９積分点
\begin{align} S'=&w_1 \left\{w_1 f(x(u_1),y(v_1)) + w_2 f(x(u_1),y(v_2)) + w_3 f(x(u_1),y(v_3))\right\}\\ +&w_2 \left\{w_1 f(x(u_2),y(v_1)) + w_2 f(x(u_2),y(v_2)) + w_3 f(x(u_2),y(v_3))\right\}\\ +&w_3 \left\{w_1 f(x(u_3),y(v_1)) + w_2 f(x(u_3),y(v_2)) + w_3 f(x(u_3),y(v_3))\right\} \end{align} 

i  j              wi              wj               ui              vj
1  1  0.555555555556  0.555555555556 -0.774596669241 -0.774596669241
1  2  0.555555555556  0.888888888889 -0.774596669241  0.000000000000
1  3  0.555555555556  0.555555555556 -0.774596669241  0.774596669241
2  1  0.888888888889  0.555555555556  0.000000000000 -0.774596669241
2  2  0.888888888889  0.888888888889  0.000000000000  0.000000000000
2  3  0.888888888889  0.555555555556  0.000000000000  0.774596669241
3  1  0.555555555556  0.555555555556  0.774596669241 -0.774596669241
3  2  0.555555555556  0.888888888889  0.774596669241  0.000000000000
3  3  0.555555555556  0.555555555556  0.774596669241  0.774596669241


inserted by FC2 system