1! 内部手続き
2module mod_integrator
3 implicit none
4
5contains
6 !
7 ! 台形公式による数値積分
8 !
9 function trapezoid(f, dx) result(ret)
10 implicit none
11 real(8), intent(in) :: f(:)
12 real(8), intent(in) :: dx
13 real(8) :: ret
14
15 integer :: i, n
16
17 n = size(f)
18
19 ret = 0.5_8 * (f(1) + f(n))
20 do i = 2, n - 1
21 ret = ret + f(i)
22 enddo
23 ret = ret * dx
24
25 endfunction trapezoid
26
27 !
28 ! Simpson公式による数値積分
29 !
30 function simpson(f, dx) result(ret)
31 implicit none
32 real(8), intent(in) :: f(:)
33 real(8), intent(in) :: dx
34 real(8) :: ret
35
36 integer :: i, n
37
38 n = size(f)
39
40 ! 端点を含めた配列サイズが奇数でなければエラー
41 if(mod(n, 2) == 0) then
42 write(*, *) 'array size must be odd'
43 stop
44 endif
45
46 ret = f(1) + f(n)
47 ! even
48 do i = 2, n - 1, 2
49 ret = ret + 4.0_8 * f(i)
50 enddo
51 ! odd
52 do i = 3, n - 2, 2
53 ret = ret + 2.0_8 * f(i)
54 enddo
55 ret = ret * dx / 3.0_8
56
57 endfunction simpson
58endmodule mod_integrator
59
60program sample
61 use mod_integrator ! モジュールを参照
62 implicit none
63
64 integer :: i, nmax
65 real(8) :: x1, x2, dx, integral, pi
66 real(8), allocatable :: y(:)
67
68 ! 真の値
69 pi = atan(1.0_8) * 4.0_8
70
71 write(*, fmt='(a)', advance='no') 'Input number of segments : '
72 read(*, *) nmax
73
74 if(mod(nmax, 2) /= 0) then
75 write(*, *) 'Number of segments must be an even number'
76 stop
77 endif
78
79 x1 = 0.0_8
80 x2 = 1.0_8
81 dx = (x2 - x1) / nmax
82
83 ! 配列サイズはnmax+1
84 allocate(y(nmax + 1))
85
86 ! 関数を評価して配列に代入
87 do i = 1, nmax + 1
88 y(i) = f(x1 + dx * real(i - 1, 8))
89 enddo
90
91 ! 台形公式による数値積分
92 integral = trapezoid(y, dx)
93 write(*, fmt='(a16, " : ", f18.15, " (rel. error = ", e18.12, ")")') &
94 & "Trapezoidal rule", integral, abs(integral / pi - 1.0_8)
95
96 ! Simpsonの公式による数値積分
97 integral = simpson(y, dx)
98 write(*, fmt='(a16, " : ", f18.15, " (rel. error = ", e18.12, ")")') &
99 & "Simpson's rule", integral, abs(integral / pi - 1.0_8)
100
101 deallocate(y)
102
103 stop
104contains
105 ! 被積分関数
106 function f(x) result(ret)
107 implicit none
108 real(8), intent(in) :: x
109 real(8) :: ret
110
111 ret = 4.0_8 / (1.0_8 + x**2)
112
113 return
114 endfunction f
115endprogram sample