1! 二分法のモジュール
2module mod_bisection
3 implicit none
4 private
5
6 integer, parameter :: default_maxit = 50
7 real(8), parameter :: default_tol = 1.0e-8_8
8
9 public :: bisection
10
11contains
12 ! 二分法により与えられた方程式の解を求める
13 subroutine bisection(f, x1, x2, error, status, maxit, tol)
14 implicit none
15 real(8), intent(inout) :: x1, x2
16 real(8), intent(out) :: error
17 integer, intent(out) :: status
18 integer, intent(in), optional :: maxit
19 real(8), intent(in), optional :: tol
20 ! 引数として関数を受け取る
21 interface
22 function f(x) result(y)
23 real(8), intent(in) :: x
24 real(8) :: y
25 endfunction f
26 endinterface
27
28 integer :: i, n
29 real(8) :: x, y, sig, tolerance
30
31 ! 最大の反復回数
32 if(.not. present(maxit)) then
33 n = default_maxit
34 else
35 n = maxit
36 endif
37
38 ! 許容誤差
39 if(.not. present(tol)) then
40 tolerance = default_tol
41 else
42 tolerance = tol
43 endif
44
45 ! x1 < x2 とする
46 if(x1 >= x2) then
47 x = x1
48 x1 = x2
49 x2 = x1
50 endif
51
52 ! f(x1)とf(x2)の符号は逆
53 if(f(x1) * f(x2) >= 0.0) then
54 status = -1
55 return
56 endif
57
58 status = 1
59 sig = sign(1.0_8, f(x2) - f(x1))
60 do i = 1, n
61 x = (x1 + x2) * 0.5_8
62 error = abs(x - x1)
63
64 ! 収束判定
65 if(error < tolerance) then
66 status = 0
67 exit
68 endif
69
70 ! 範囲を縮小
71 y = f(x)
72 if(y * sig < 0.0) then
73 x1 = x
74 else
75 x2 = x
76 endif
77 enddo
78
79 ! 最終的な根の近似値
80 x1 = x
81 return
82 endsubroutine bisection
83
84endmodule mod_bisection