積分プログラム(長方形積分) [fortran]

こんにちは、管理人のさくらえびです。
今回の記事は、積分プログラムの初歩。長方形積分について書きます。
プログラムで面積を計算する、誰もが一番最初に通るのと思います。
長方形プログラムでは、斜めの関数や斜めの直線を、長方形の塊として捉えて計算します。
式で書くとこんな感じです。
\[
\int_a^b f(x)dx → \sum_{i=1}^{n-1} f(x_i)dx
\]
となります。ここで\(n\)はメッシュ(分割数)です。
簡単な絵で描くと、こんな感じです。

総和の形になったので、これで計算できます。
例を示します。最初にこんな感じの関数を用意します。
\[
f(x) = -x+1
\]
0から1の範囲で積分するとき、答えはを簡単に0.5と求まります。
ではこれを、長方形積分したときに、どうなるでしょうか。
下記は積分のソースコードです。

subroutine SQ_INTEG(func,a,b,mesh,s)
implicit none
double precision func,a,b,s,s0,dx,pos
integer*4 mesh,i
s0 = 0.d0
pos = a
dx = (b-a)/dble(mesh)
do i=1,mesh
s0 = s0 + dx*func(pos)
pos = pos+dx
enddo
s = s0
end subroutine

さて、これを\(f(x) = -x+1\)を0から1までの積分でテストしてみましょう。
テスト用のコードはこんな感じ

program test
implicit none
external fx
double precision x0,x1,fx,s
integer*4 mesh
x0 = 0
x1 = 1
mesh = 10
call SQ_INTEG(fx,x0,x1,mesh,s)
print *, "分割数,面積,計算誤差=",mesh,s,s-0.5d0
end program
function fx(x)
implicit none
double precision fx,x
fx = -x + 1
end function

結果はこんな感じ

誤差は、0.05になっていますね。10(分割数)×0.1(分割したときの微小x)×0.1(\(f(0.1)-f(0)\)の差)/2=0.05が乗ったことになります。
大きいですね。じゃあ、メッシュを９回くらい10倍ずつしてみた結果を見てみましょう。
なぜ9回かというと、整数型の最大値を超えてしまうからです。
下記は結果です。

数が増えるごとに、誤差は一桁ずつ小さくなってますね。
ただ、整数値にも限りがあるので、整数値の最大までとなります。
倍精度実数を使えば、もっといけますが、倍精度実数の桁も誤差の許容範囲があるので、\(10^{-13}\)位が限界ですね。
それだけやっても、結局計算時間ばかり使うので、あまりうれしいことはありません。
長方形積分は、演算時間と精度がトレードオフになることが確認できました
by さくらえび