1program sample
2 implicit none
3
4 integer, parameter :: nmax = 50
5 real(8), parameter :: tolerance = 1.0e-6
6
7 logical :: status
8 integer :: n
9 real(8) :: x, y, x1, x2, sig
10
11 write(*, fmt='(a)', advance='no') 'Input two initial guesses : '
12 read(*, *) x1, x2
13
14 ! x1 < x2 となるように交換
15 if(x1 >= x2) then
16 x = x1
17 x1 = x2
18 x2 = x
19 endif
20
21 ! f(x1)とf(x2)の符号は逆でなければならない
22 if(f(x1) * f(x2) >= 0.0) then
23 write(*, *) 'Signs of two function values must be opposite'
24 stop
25 endif
26
27 sig = sign(1.0_8, f(x2) - f(x1))
28 status = .false.
29 do n = 1, nmax
30 x = (x1 + x2) * 0.5_8
31 y = f(x)
32
33 ! 収束判定
34 if(abs(x2 - x1) < tolerance) then
35 status = .true.
36 write(*, fmt='("Converged at step ", i4)') n
37 exit
38 else
39 write(*, fmt='("Error at step ", i4, " = ", e20.8)') n, abs(y)
40 endif
41
42 ! 次の値を推定
43 if(y * sig < 0.0) then
44 x1 = x
45 else
46 x2 = x
47 endif
48 enddo
49
50 ! 結果の表示
51 if(status .eqv. .false.) then
52 write(*, fmt='("Failed to find a root after ", i3, " iterations")') nmax
53 endif
54
55 write(*, fmt='("Final approximation = ", e20.8)') x
56 write(*, fmt='("Final error = ", e20.8)') abs(y)
57
58 stop
59contains
60 ! 方程式
61 function f(x) result(ret)
62 implicit none
63 real(8), intent(in) :: x
64 real(8) :: ret
65
66 ret = x - cos(x)
67
68 return
69 endfunction f
70
71endprogram sample