1program answer
2 implicit none
3
4 integer :: i, n
5 real(8) :: x1, x2, integral, err1, err2
6
7 ! 真の値
8 integral = 1.0_8
9
10 x1 = 0.0_8
11 x2 = 1.0_8
12
13 do i = 1, 16
14 n = 2**i
15 err1 = abs(trapezoid(n, x1, x2) / integral - 1.0_8)
16 err2 = abs(simpson(n, x1, x2) / integral - 1.0_8)
17 write(*, fmt='(i8, 1x, e12.4, 1x, e12.4)') n, err1, err2
18 enddo
19
20 stop
21contains
22 ! f(x)
23 function f(x) result(ret)
24 implicit none
25 real(8), intent(in) :: x
26 real(8) :: ret
27
28 real(8), parameter :: factor = 1.0_8 / atan(1.0_8)
29
30 ret = factor / (1.0_8 + x**2)
31 !ret = 2*x
32 !ret = 3*x**2
33 !ret = 4*x**3
34 !ret = 5*x**4
35 !ret = 6*x**5
36
37 return
38 endfunction f
39
40 ! 台形公式
41 function trapezoid(n, a, b) result(ret)
42 implicit none
43 integer, intent(in) :: n
44 real(8), intent(in) :: a, b
45 real(8) :: ret
46
47 integer :: i
48 real(8) :: h
49
50 h = (b - a) / n
51
52 ret = 0.5_8 * (f(a) + f(b))
53 do i = 1, n - 1
54 ret = ret + f(a + h * real(i, 8))
55 enddo
56
57 ret = ret * h
58
59 endfunction trapezoid
60
61 ! Simpsonの公式
62 function simpson(n, a, b) result(ret)
63 implicit none
64 integer, intent(in) :: n
65 real(8), intent(in) :: a, b
66 real(8) :: ret
67
68 integer :: i
69 real(8) :: h
70
71 h = (b - a) / n
72
73 ret = f(x1) + f(x2)
74 ! odd
75 do i = 1, n - 1, 2
76 ret = ret + 4.0_8 * f(a + h * real(i, 8))
77 enddo
78 ! even
79 do i = 2, n - 2, 2
80 ret = ret + 2.0_8 * f(a + h * real(i, 8))
81 enddo
82
83 ret = ret * h / 3.0_8
84
85 endfunction simpson
86
87endprogram answer