Fortran - 数値積分(台形則/シンプソン則による定積分)!

Updated:


Fortran 95 で、数値積分(台形則/シンプソン則による定積分)行ってみました。

0. 前提条件

  • LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
  • GCC 6.3.0 (GFortran 6.3.0) でのコンパイルを想定。

1. 数値積分(台形則/シンプソン則による定積分)とは?

当ブログ過去記事を参照のこと。

2. 想定する被積分関数

想定する被積分関数は \(f(x) = \sqrt{4 - x^{2}}\ \ (但し、 0.0 \leq x \leq 2.0\))

DEFINITE_INTEGRAL

3. ソースコードの作成

  • 台形則/シンプソン則計算部分はモジュール化している。
  • 小さなファイルなので、モジュール部分を別ファイルには分けていない。

File: definite_integral.f95

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
!****************************************************
! 定積分(台形則、シンプソン則)
! * f = sqrt(4.0 - x**2)  (0.0 <= x <= 2.0)
!
! date          name            version
! 2018.10.13    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module integrator
  implicit none
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
  real(DP),    parameter :: PI = atan(1.0_8) * 4.0_DP  ! 円周率

contains
  ! 数値積分(台形則)
  !
  ! :param  real(8) f(:)
  ! :param  real(8)   dx
  ! :return real(8)  ret
  function trapezoid(f, dx) result(ret)
    implicit none
    real(DP), intent(in) :: f(:)
    real(DP), intent(in) :: dx
    real(DP)    :: ret
    integer(SP) :: i, n

    n = size(f)
    ret = 0.5_DP * (f(1) + f(n))
    do i = 2, n - 1
       ret = ret + f(i)
    end do
    ret = ret * dx
  end function trapezoid

  ! 数値積分(シンプソン則)
  !
  ! :param  real(8)  f(:)
  ! :param  real(8)    dx
  ! :return real(8)   ret
  function simpson(f, dx) result(ret)
    implicit none
    real(DP), intent(in) :: f(:)
    real(DP), intent(in) :: dx
    real(DP)    :: ret
    integer(SP) :: i, n

    n = size(f)
    ! 端点を含めた配列サイズが奇数でなければ終了
    if (mod(n, 2) == 0) then
       write(*,*) '配列サイズが奇数でない!'
       stop
    end if
    ret = f(1) + f(n)
    ! (偶数)
    do i = 2, n - 1, 2
       ret = ret + 4.0_8 * f(i)
    end do
    ! (奇数)
    do i = 3, n - 2, 2
       ret = ret + 2.0_DP * f(i)
    end do
    ret = ret * dx / 3.0_DP
  end function simpson
end module integrator

program definite_integral
  use integrator
  implicit none
  integer(SP) :: n, nmax
  real(DP)    :: x_1, x_2, dx, integral
  real(DP), allocatable :: f_1(:), f_2(:)

  write(*, fmt='(A)', advance='no') '区間の数 : '
  read(*,*) nmax

  ! 配列用メモリ確保
  allocate(f_1(0:nmax))
  allocate(f_2(0:nmax))

  x_1 = 0.0_DP
  x_2 = 2.0_DP
  dx = (x_2 - x_1) / nmax

  ! 数値積分(台形則)
  do n = 0, nmax
    f_1(n) = f(x_1 + dx * real(n, 8))
  end do
  integral = trapezoid(f_1, dx)
  write (*, fmt='(a15, ": ", f18.15, " (誤差= ", e18.12, ")")') &
    & "台形則", integral, abs(integral / PI - 1.0_DP)

  ! 数値積分(シンプソン則)
  do n = 0, nmax
    f_2(n) = f(x_1 + dx * real(n, 8))
  end do
  integral = simpson(f_2, dx)
  write(*, fmt='(a18, ": ", f18.15, " (誤差= ", e18.12, ")")') &
       & "シンプソン則", integral, abs(integral / PI - 1.0_DP)

  ! 配列用メモリ解放
  deallocate(f_1)
  deallocate(f_2)

  stop
contains
  ! 被積分関数
  real(DP) function f(x)
    implicit none
    real(DP), intent(in) :: x

    f = sqrt(4.0_DP - x * x)
  end function f
end program definite_integral

4. ソースコードのコンパイル

$ gfortran -o definite_integral definite_integral.f95

5. 動作確認

区間を分割する数を問われるので、入力する。(シンプソン則にも対応するために偶数で指定する)

$ ./definite_integral
区間の数 : 100
      台形則:  3.140417031779045 (誤差= 0.374212044774E-03)
シンプソン則:  3.141133205339227 (誤差= 0.146246920345E-03)

以上、





 

Sponsored Link

 

Comments