Simpson法
比較的小さい積分領域 $x_a \leqq x \leqq x_b$ を領域を考えたとき,一次元のSimpson法は以下のとおり表現できた.
\begin{equation}
\int_{x_a}^{x_b} f(x) dx = \cfrac{x_b-x_a}{6} \left\{f(x_a) + 4 f \left(\cfrac{x_a+x_b}{2}\right) + f(x_b)\right\}
\end{equation}
上式を二次元に拡張することを考える.積分領域は $x_a \leqq x \leqq x_b$,$y_a \leqq y \leqq y_b$ の長方形領域を考え, 表現の単純化のため,積分領域の中間点 $x_c=(x_a+x_b)/2$,$y_c=(y_a+y_b)/2$ を定義しておく.
まず $y$ 方向の3測線を考え,それぞれの積分値を $s_1$,$s_2$,$s_3$ とし,これを $x$ 方向で積分することを考えると,以下のように2重積分値 $S$ が定められる
\begin{align}
&s_1=\cfrac{y_b-y_a}{6}\left\{f(x_a,y_a)+4 f(x_a,y_c)+f(x_a,y_b)\right\}\\
&s_2=\cfrac{y_b-y_a}{6}\left\{f(x_c,y_a)+4 f(x_c,y_c)+f(x_c,y_b)\right\}\\
&s_3=\cfrac{y_b-y_a}{6}\left\{f(x_b,y_a)+4 f(x_b,y_c)+f(x_b,y_b)\right\}\\
&S=\cfrac{x_b-x_a}{6}(s_1 + 4 s_2 + s_3)
\end{align}
Gauss-Legendre積分
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\int 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\}$ と置き,Simpson法との整合を考え,9積分点モデルとすれば,以下のように $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)) + 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 | ai | bj |
---|---|---|---|---|---|
1 | 1 | 0.555555555555555 | 0.555555555555555 | -0.774596669241483 | -0.774596669241483 |
1 | 2 | 0.555555555555555 | 0.888888888888889 | -0.774596669241483 | 0.000000000000000 |
1 | 3 | 0.555555555555555 | 0.555555555555555 | -0.774596669241483 | 0.774596669241483 |
2 | 1 | 0.888888888888889 | 0.555555555555555 | 0.000000000000000 | -0.774596669241483 |
2 | 2 | 0.888888888888889 | 0.888888888888889 | 0.000000000000000 | 0.000000000000000 |
2 | 3 | 0.888888888888889 | 0.555555555555555 | 0.000000000000000 | 0.774596669241483 |
3 | 1 | 0.555555555555555 | 0.555555555555555 | 0.774596669241483 | -0.774596669241483 |
3 | 2 | 0.555555555555555 | 0.888888888888889 | 0.774596669241483 | 0.000000000000000 |
3 | 3 | 0.555555555555555 | 0.555555555555555 | 0.774596669241483 | 0.774596669241483 |
Programs
Programs
Filename | Description |
---|---|
f90_DSIMP.txt | program for 2-dim Simpson integration |
f90_DGINT3.txt | program for 2-dim. Gauss-Legendre integration |
a_GMT.txt | Script for GMT drawing |
解の収束状況作図事例
以下に計算事例を示す.積分した式は以下のとおり.
\begin{equation}
\int_0^1 \int_0^1 \exp(x+y) dx dy = (e-1)^2
\end{equation}
Simpson | |||
---|---|---|---|
nx | ny | Estimated value | Error |
1 | 1 | 0.295448365943E+01 | +0.199121741797E-02 |
2 | 2 | 0.295261964250E+01 | +0.127200490735E-03 |
3 | 3 | 0.295251767147E+01 | +0.252294546619E-04 |
4 | 4 | 0.295250043629E+01 | +0.799428018006E-05 |
5 | 5 | 0.295249571866E+01 | +0.327664674993E-05 |
7 | 7 | 0.295249329545E+01 | +0.853435002668E-06 |
10 | 10 | 0.295249264699E+01 | +0.204973195750E-06 |
20 | 20 | 0.295249245483E+01 | +0.128136812272E-07 |
30 | 30 | 0.295249244454E+01 | +0.253120102656E-08 |
40 | 40 | 0.295249244281E+01 | +0.800902455467E-09 |
50 | 50 | 0.295249244234E+01 | +0.328065574706E-09 |
70 | 70 | 0.295249244210E+01 | +0.853996873218E-10 |
100 | 100 | 0.295249244203E+01 | +0.204396499726E-10 |
200 | 200 | 0.295249244201E+01 | +0.113864473406E-11 |
300 | 300 | 0.295249244201E+01 | +0.287325718773E-12 |
400 | 400 | 0.295249244201E+01 | -0.339284156325E-12 |
500 | 500 | 0.295249244201E+01 | +0.329514193709E-12 |
700 | 700 | 0.295249244201E+01 | +0.198507876803E-12 |
1000 | 1000 | 0.295249244201E+01 | +0.244249065418E-13 |
Gauss-Legendre | |||
---|---|---|---|
nx | ny | Estimated value | Error |
1 | 1 | 0.295248960999E+01 | -0.283202512019E-05 |
2 | 2 | 0.295249239663E+01 | -0.453798536526E-07 |
3 | 3 | 0.295249243801E+01 | -0.400277144763E-08 |
4 | 4 | 0.295249244130E+01 | -0.713587411383E-09 |
5 | 5 | 0.295249244183E+01 | -0.187207582769E-09 |
7 | 7 | 0.295249244199E+01 | -0.248818743387E-10 |
10 | 10 | 0.295249244201E+01 | -0.293143287422E-11 |
20 | 20 | 0.295249244201E+01 | -0.492939022934E-13 |
30 | 30 | 0.295249244201E+01 | +0.444089209850E-15 |
40 | 40 | 0.295249244201E+01 | +0.000000000000E+00 |
50 | 50 | 0.295249244201E+01 | +0.399680288865E-14 |
70 | 70 | 0.295249244201E+01 | -0.164313007645E-13 |
100 | 100 | 0.295249244201E+01 | -0.310862446895E-13 |
200 | 200 | 0.295249244201E+01 | -0.123456800338E-12 |
300 | 300 | 0.295249244201E+01 | +0.568434188608E-13 |
400 | 400 | 0.295249244201E+01 | -0.398792110445E-12 |
500 | 500 | 0.295249244201E+01 | +0.300648395068E-12 |
700 | 700 | 0.295249244201E+01 | +0.193178806285E-12 |
1000 | 1000 | 0.295249244201E+01 | +0.244249065418E-13 |