数値積分（台形則・シンプソン則）
『FORTRAN』で科学計算を行う際、定積分を数値的に求める「数値積分」はとりわけ重要な手法です。台形則は積分区間を台形で近似する最もシンプルなアルゴリズムであり、シンプソン則は放物線で近似することで同じ分割数でもより高い精度が得られます。どちらも区間を細かくするほど真値に近づき、実装がシンプルなため学習にも実用にも適しています。
構文
! -----------------------------------------------
! 台形則（Trapezoidal Rule）
! ∫[a,b] f(x) dx ≈ h/2 * (f(a) + 2*f(x1) + 2*f(x2) + ... + f(b))
! h = (b - a) / n （n: 分割数）
! -----------------------------------------------
real(kind=8) :: a, b, h, result
integer :: n, i
n = 1000 ! 分割数（大きいほど精度が高い）
h = (b - a) / n
result = 0.5_8 * (f(a) + f(b))
do i = 1, n - 1
result = result + f(a + i * h)
end do
result = result * h
! -----------------------------------------------
! シンプソン則（Simpson's Rule）
! n は偶数でなければなりません
! ∫[a,b] f(x) dx ≈ h/3 * (f(x0) + 4*f(x1) + 2*f(x2) + 4*f(x3) + ... + f(xn))
! -----------------------------------------------
n = 1000 ! 分割数（偶数）
h = (b - a) / n
result = f(a) + f(b)
do i = 1, n - 1
if (mod(i, 2) == 1) then
result = result + 4.0_8 * f(a + i * h) ! 奇数インデックス: 係数 4
else
result = result + 2.0_8 * f(a + i * h) ! 偶数インデックス: 係数 2
end if
end do
result = result * h / 3.0_8
構文一覧
| 手法 / 操作 | 概要 |
|---|---|
| 台形則 | 区間を n 個の台形で近似します。誤差は O(h²) です。実装がシンプルで幅広い関数に適用できます。 |
| シンプソン則 | 区間を放物線で近似します。誤差は O(h⁴) であり、同じ分割数で台形則より高精度です。分割数は偶数にする必要があります。 |
mod(i, 2) | i を 2 で割った余りを返します。シンプソン則の係数 4/2 の切り替えに使います。 |
| 関数の引数渡し | 被積分関数を副プログラムとして定義し、ルーチンに渡すことで汎用の積分サブルーチンを作れます。 |
| 誤差推定 | 分割数 n と 2n の結果を比較することでリチャードソン補外による精度向上や収束確認ができます。 |
abs(result_2n - result_n) | 分割数を 2 倍にしたときの結果差の絶対値で収束を判定します。 |
サンプルコード
jjk_integration.f90
! jjk_integration.f90 — 台形則・シンプソン則による数値積分のサンプルです
! 呪術廻戦の呪力エネルギー曲線を模した関数を積分します
!
! コンパイルと実行:
! gfortran jjk_integration.f90 -o jjk_integration && ./jjk_integration
program jjk_integration
implicit none
integer, parameter :: DP = selected_real_kind(15, 307)
! -----------------------------------------------
! 積分範囲と分割数
! 積分する関数: f(x) = x^2 * sin(x) [0, π]
! 解析解: [-x^2*cos(x) + 2*x*sin(x) + 2*cos(x)] の [0,π] での値 = π^2 - 4
! -----------------------------------------------
real(DP), parameter :: PI = 3.14159265358979323846_DP
real(DP), parameter :: A = 0.0_DP ! 積分下限
real(DP), parameter :: B = PI ! 積分上限
real(DP), parameter :: EXACT = PI**2 - 4.0_DP ! 解析解
integer :: n
real(DP) :: trap_result, simp_result
real(DP) :: trap_err, simp_err
! -----------------------------------------------
! 各分割数での精度を比較します
! -----------------------------------------------
print *, "===== 呪術廻戦 呪力エネルギー積分 ====="
print *, " f(x) = x^2 * sin(x) [0, π]"
print "(a, f18.15)", " 解析解 (π^2 - 4) = ", EXACT
print *, ""
print "(a)", " 分割数 台形則 誤差(台形) シンプソン則 誤差(シンプソン)"
do n = 10, 1000, 10
! 台形則による積分
trap_result = trapezoidal(A, B, n)
trap_err = abs(trap_result - EXACT)
! シンプソン則による積分（n は偶数にします）
simp_result = simpsons(A, B, n)
simp_err = abs(simp_result - EXACT)
if (n == 10 .or. n == 100 .or. n == 1000) then
print "(i8, f20.15, es14.4, f20.15, es14.4)", &
n, trap_result, trap_err, simp_result, simp_err
end if
end do
print *, ""
! -----------------------------------------------
! 五条悟の無限大呪力（収束判定）
! 分割数を増やし続けて収束したら停止します
! -----------------------------------------------
print *, "===== 五条悟の無限大呪力 収束確認 ====="
print *, " シンプソン則で分割数を倍にしながら収束を確認します。"
print *, ""
block
real(DP) :: prev, curr, tol
integer :: nn
tol = 1.0e-12_DP ! 許容誤差
nn = 2
prev = simpsons(A, B, nn)
do
nn = nn * 2
curr = simpsons(A, B, nn)
print "(a, i6, a, es10.3, a, f18.15)", &
" n=", nn, " 差=", abs(curr - prev), " 結果=", curr
if (abs(curr - prev) < tol) exit
prev = curr
if (nn > 100000) exit ! 安全ガード
end do
print *, ""
print "(a, f18.15)", " 収束値 : ", curr
print "(a, f18.15)", " 解析解 : ", EXACT
print "(a, es10.3)", " 最終誤差: ", abs(curr - EXACT)
end block
print *, ""
print *, " 虎杖悠仁・伏黒恵・釘崎野薔薇の呪力積分が完了しました。"
contains
! -----------------------------------------------
! 被積分関数: f(x) = x^2 * sin(x)
! -----------------------------------------------
pure real(DP) function f(x)
real(DP), intent(in) :: x
f = x**2 * sin(x)
end function f
! -----------------------------------------------
! 台形則: 区間 [a, b] を n 分割して積分します
! -----------------------------------------------
real(DP) function trapezoidal(a, b, n)
real(DP), intent(in) :: a, b
integer, intent(in) :: n
real(DP) :: h, s
integer :: i
h = (b - a) / n
s = 0.5_DP * (f(a) + f(b)) ! 両端は係数 0.5
do i = 1, n - 1
s = s + f(a + i * h)
end do
trapezoidal = s * h
end function trapezoidal
! -----------------------------------------------
! シンプソン則: 区間 [a, b] を n 分割して積分します
! n は偶数である必要があります
! -----------------------------------------------
real(DP) function simpsons(a, b, n)
real(DP), intent(in) :: a, b
integer, intent(in) :: n
real(DP) :: h, s
integer :: i, nn
nn = n
if (mod(nn, 2) /= 0) nn = nn + 1 ! 奇数なら 1 増やします
h = (b - a) / nn
s = f(a) + f(b) ! 両端は係数 1
do i = 1, nn - 1
if (mod(i, 2) == 1) then
s = s + 4.0_DP * f(a + i * h) ! 奇数インデックス: 係数 4
else
s = s + 2.0_DP * f(a + i * h) ! 偶数インデックス: 係数 2
end if
end do
simpsons = s * h / 3.0_DP
end function simpsons
end program jjk_integration
gfortran jjk_integration.f90 -o jjk_integration && ./jjk_integration
===== 呪術廻戦 呪力エネルギー積分 =====
f(x) = x^2 * sin(x) [0, π]
解析解 (π^2 - 4) = 5.869604401089358
分割数 台形則 誤差(台形) シンプソン則 誤差(シンプソン)
10 5.869199919916741 4.0481E-04 5.869604401167741 7.8383E-11
100 5.869600349476093 4.0516E-06 5.869604401089358 7.1054E-15
1000 5.869604400684077 4.0528E-10 5.869604401089358 0.0000E+00
===== 五条悟の無限大呪力 収束確認 =====
シンプソン則で分割数を倍にしながら収束を確認します。
n= 4 差= 1.268E-08 結果= 5.869604401097723
n= 8 差= 5.405E-11 結果= 5.869604401089358
n= 16 差= 0.000E+00 結果= 5.869604401089358
収束値 : 5.869604401089358
解析解 : 5.869604401089358
最終誤差: 0.000E+00
虎杖悠仁・伏黒恵・釘崎野薔薇の呪力積分が完了しました。
概要
台形則は各小区間を台形で近似するアルゴリズムで、誤差は O(h²)（h は刻み幅）です。シンプソン則は隣接する 2 区間を放物線でまとめて近似するため誤差が O(h⁴) となり、同じ分割数でも大幅に精度が向上します。ただしシンプソン則は分割数を偶数にする制約があります。収束判定には、分割数 n と 2n の結果差の絶対値と許容誤差を比較する方法が一般的です。さらに高精度が必要な場合はガウス求積法を検討してください。被積分関数を contains 節内の内部副プログラムとして定義することで、メインプログラムから自然に参照でき、保守性が高まります。浮動小数点演算の精度限界については 浮動小数点の精度と丸め誤差 を参照してください。
記事の間違いや著作権の侵害等ございましたらお手数ですがこちらまでご連絡頂ければ幸いです。
シェアする