1program sample
2 implicit none
3
4 integer :: n, nmax
5 real(8) :: x1, x2, dx, f1, f2, integral, pi
6
7 ! 真の値
8 pi = atan(1.0_8) * 4.0_8
9
10 write(*, fmt='(a)', advance='no') 'Input number of segments : '
11 read(*, *) nmax
12
13 if(mod(nmax, 2) /= 0) then
14 write(*, *) 'Number of segments must be an even number'
15 stop
16 endif
17
18 x1 = 0.0_8
19 x2 = 1.0_8
20 dx = (x2 - x1) / nmax
21
22 !
23 ! 台形公式による数値積分
24 !
25 integral = 0.5_8 * (f(x1) + f(x2))
26 do n = 1, nmax - 1
27 integral = integral + f(x1 + dx * real(n, 8))
28 enddo
29 integral = integral * dx
30
31 write(*, fmt='(a16, " : ", f18.15, " (rel. error = ", e18.12, ")")') &
32 & "Trapezoidal rule", integral, abs(integral / pi - 1.0_8)
33
34 !
35 ! Simpsonの公式による数値積分
36 !
37 integral = f(x1) + f(x2)
38 ! odd
39 do n = 1, nmax - 1, 2
40 integral = integral + 4.0_8 * f(x1 + dx * real(n, 8))
41 enddo
42 ! even
43 do n = 2, nmax - 2, 2
44 integral = integral + 2.0_8 * f(x1 + dx * real(n, 8))
45 enddo
46 integral = integral * dx / 3.0_8
47
48 write(*, fmt='(a16, " : ", f18.15, " (rel. error = ", e18.12, ")")') &
49 & "Simpson's rule", integral, abs(integral / pi - 1.0_8)
50
51 stop
52contains
53 ! 被積分関数
54 function f(x) result(ret)
55 implicit none
56 real(8), intent(in) :: x
57 real(8) :: ret
58
59 ret = 4.0_8 / (1.0_8 + x**2)
60
61 return
62 endfunction f
63endprogram sample