1! Newton法のモジュール
2module mod_newton
3 implicit none
4 private
5
6 integer, parameter :: default_maxit = 50
7 real(8), parameter :: default_tol = 1.0e-8_8
8 real(8), parameter :: epsilon = 1.0e-308_8
9
10 ! オーバーロード
11 interface newton
12 module procedure newton_func, newton_sub
13 endinterface newton
14
15 public :: newton
16
17contains
18 ! Newton法 (fおよびdf/dxをそれぞれ返す関数を引数として受け取る)
19 subroutine newton_func(f, df, x, error, status, maxit, tol)
20 implicit none
21 real(8), intent(inout) :: x
22 real(8), intent(out) :: error
23 integer, intent(out) :: status
24 integer, intent(in), optional :: maxit
25 real(8), intent(in), optional :: tol
26 ! 引数としてサブルーチンを受け取る
27 interface
28 ! 関数値
29 real(8) function f(x)
30 real(8), intent(in) :: x
31 endfunction f
32 ! 微分値
33 real(8) function df(x)
34 real(8), intent(in) :: x
35 endfunction df
36 endinterface
37
38 integer :: i, n
39 real(8) :: dx, y, dy, tolerance
40
41 ! 最大の反復回数
42 if(.not. present(maxit)) then
43 n = default_maxit
44 else
45 n = maxit
46 endif
47
48 ! 許容誤差
49 if(.not. present(tol)) then
50 tolerance = default_tol
51 else
52 tolerance = tol
53 endif
54
55 status = 1
56 do i = 1, n
57 y = f(x)
58 dy = df(x)
59 dx = -y / (dy + epsilon)
60 x = x + dx
61 error = abs(dx)
62
63 ! 収束判定
64 if(error < tolerance) then
65 status = 0
66 exit
67 endif
68 enddo
69
70 return
71 endsubroutine newton_func
72
73 ! Newton法 (fおよびdf/dxを返すサブルーチンを引数として受け取る)
74 subroutine newton_sub(fdf, x, error, status, maxit, tol)
75 implicit none
76 real(8), intent(inout) :: x
77 real(8), intent(out) :: error
78 integer, intent(out) :: status
79 integer, intent(in), optional :: maxit
80 real(8), intent(in), optional :: tol
81 ! 引数としてサブルーチンを受け取る
82 interface
83 subroutine fdf(x, f, df)
84 real(8), intent(in) :: x
85 real(8), intent(out) :: f, df
86 endsubroutine fdf
87 endinterface
88
89 integer :: i, n
90 real(8) :: dx, y, dy, tolerance
91
92 ! 最大の反復回数
93 if(.not. present(maxit)) then
94 n = default_maxit
95 else
96 n = maxit
97 endif
98
99 ! 許容誤差
100 if(.not. present(tol)) then
101 tolerance = default_tol
102 else
103 tolerance = tol
104 endif
105
106 status = 1
107 do i = 1, n
108 call fdf(x, y, dy)
109 dx = -y / (dy + epsilon)
110 x = x + dx
111 error = abs(dx)
112
113 ! 収束判定
114 if(error < tolerance) then
115 status = 0
116 exit
117 endif
118 enddo
119
120 return
121 endsubroutine newton_sub
122
123endmodule mod_newton
124
125module mod_equation
126 implicit none
127contains
128 ! 関数値
129 function f(x) result(y)
130 real(8), intent(in) :: x
131 real(8) :: y
132
133 y = x - cos(x)
134 endfunction f
135
136 ! 微分値
137 function df(x) result(y)
138 real(8), intent(in) :: x
139 real(8) :: y
140
141 y = 1 + sin(x)
142 endfunction df
143
144 ! サブルーチンで2つの値を同時に返す
145 subroutine fdf(x, y, dy)
146 implicit none
147 real(8), intent(in) :: x
148 real(8), intent(out) :: y, dy
149
150 y = f(x)
151 dy = df(x)
152
153 return
154 endsubroutine fdf
155endmodule mod_equation
156
157program answer
158 use mod_newton
159 use mod_equation
160 implicit none
161
162 integer :: status
163 real(8) :: x, err
164
165 x = 1.0
166
167 ! 関数f, dfを引数として渡す
168 call newton(f, df, x, err, status, 20, 1.0e-6_8)
169
170 ! サブルーチンfdfを引数として渡す
171 !call newton(fdf, x, err, status, 20, 1.0e-6_8)
172
173 select case(status)
174 case(0)
175 write(*, fmt='(a)') 'Iteration converged !'
176 case(1)
177 write(*, fmt='(a)') 'Iteration did not converge !'
178 case(-1)
179 write(*, fmt='(a)') 'Error'
180 stop
181 case default
182 write(*, fmt='(a)') 'Unknown error'
183 stop
184 endselect
185
186 write(*, fmt='("Solution = ", e18.12, ", Error = ", e8.2)') x, err
187
188 stop
189endprogram answer