1program answer
2 implicit none
3
4 integer, parameter :: nmax = 20
5 real(8), parameter :: epsilon = 1.0e-308_8
6 real(8), parameter :: rm = sqrt(1836.0_8)
7
8 integer :: i
9 real(8) :: x1, x2
10 real(8) :: err1(nmax), err2(nmax), err3(nmax)
11
12 ! bisection method
13 x1 = -2.6_8
14 x2 = -2.4_8
15 call bisection(x1, x2, err1)
16
17 ! secant method
18 x1 = 1.0_8
19 x2 = 0.0_8
20 call secant(x1, x2, err2)
21
22 ! newton method
23 x1 = 0.0_8
24 call newton(x1, err3)
25
26 do i = 1, nmax
27 write(*, '(i3, 1x, e12.4, 1x, e12.4, 1x, e12.4)') &
28 & i, err1(i), err2(i), err3(i)
29 enddo
30
31contains
32
33 ! f(x)
34 function f(x) result(y)
35 implicit none
36 real(8), intent(in) :: x
37 real(8) :: y
38
39 y = (1.0_8 - x) / rm - exp(x)
40 endfunction f
41
42 ! f'(x)
43 function df(x) result(y)
44 implicit none
45 real(8), intent(in) :: x
46 real(8) :: y
47
48 y = -1.0_8 / rm - exp(x)
49 endfunction df
50
51 ! 二分法
52 subroutine bisection(x1, x2, error)
53 implicit none
54 real(8), intent(inout) :: x1, x2
55 real(8), intent(out) :: error(nmax)
56
57 integer :: n
58 real(8) :: x, y, sig
59
60 ! x1 < x2 とする
61 if(x1 >= x2) then
62 x = x1
63 x1 = x2
64 x2 = x1
65 endif
66
67 ! f(x1)とf(x2)の符号は逆
68 if(f(x1) * f(x2) >= 0.0) then
69 write(*, *) 'The signs of f(x1) and f(x2) must be opposite'
70 return
71 endif
72
73 sig = sign(1.0_8, f(x2) - f(x1))
74 do n = 1, nmax
75 x = (x1 + x2) * 0.5_8
76 error(n) = abs(x - x1)
77
78 y = f(x)
79 if(y * sig < 0.0) then
80 x1 = x
81 else
82 x2 = x
83 endif
84 enddo
85
86 x1 = x
87 return
88 endsubroutine bisection
89
90 ! Secant method (割線法)
91 subroutine secant(x1, x2, error)
92 implicit none
93 real(8), intent(inout) :: x1, x2
94 real(8), intent(out) :: error(nmax)
95
96 integer :: n
97 real(8) :: dx, dy, y1, y2
98
99 y1 = f(x1)
100 y2 = f(x2)
101 do n = 1, nmax
102 dy = (x2 - x1) / (y2 - y1 + epsilon)
103 dx = -y2 * dy
104 error(n) = abs(dx)
105 x1 = x2
106 x2 = x2 + dx
107 y1 = y2
108 y2 = f(x2)
109 enddo
110
111 x1 = x2
112 return
113 endsubroutine secant
114
115 ! Newton法
116 subroutine newton(x, error)
117 implicit none
118 real(8), intent(inout) :: x
119 real(8), intent(out) :: error(nmax)
120
121 integer :: n
122 real(8) :: dx
123
124 do n = 1, nmax
125 dx = -f(x) / (df(x) + epsilon)
126 x = x + dx
127 error(n) = abs(dx)
128 enddo
129
130 return
131 endsubroutine newton
132
133endprogram answer