rootstest.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: rootstest.f90 #
5!# #
6!# Copyright (C) 2016 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# #
9!# This program is free software; you can redistribute it and/or modify #
10!# it under the terms of the GNU General Public License as published by #
11!# the Free Software Foundation; either version 2 of the License, or (at #
12!# your option) any later version. #
13!# #
14!# This program is distributed in the hope that it will be useful, but #
15!# WITHOUT ANY WARRANTY; without even the implied warranty of #
16!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17!# NON INFRINGEMENT. See the GNU General Public License for more #
18!# details. #
19!# #
20!# You should have received a copy of the GNU General Public License #
21!# along with this program; if not, write to the Free Software #
22!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23!# #
24!#############################################################################
25
26!----------------------------------------------------------------------------!
31!----------------------------------------------------------------------------!
32PROGRAM rootstest
33 USE roots
34#include "tap.h"
35 IMPLICIT NONE
36 !--------------------------------------------------------------------------!
37 INTERFACE
38 PURE SUBROUTINE func(x,fx,plist)
39 IMPLICIT NONE
40 REAL, INTENT(IN) :: x
41 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
42 REAL, INTENT(OUT) :: fx
43 END SUBROUTINE func
44 END INTERFACE
45 INTERFACE
46 PURE SUBROUTINE funcd(x,fx,dfx,plist)
47 IMPLICIT NONE
48 REAL, INTENT(IN) :: x
49 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
50 REAL, INTENT(OUT) :: fx,dfx
51 END SUBROUTINE funcd
52 END INTERFACE
53 COMMON /test_num/ k
54 !--------------------------------------------------------------------------!
55 INTEGER, PARAMETER :: num_tests = 13
56 INTEGER, PARAMETER :: num_methods = 8
57 CHARACTER(LEN=16), PARAMETER, DIMENSION(NUM_METHODS) :: method_name = (/ &
58 " Bisection", &
59 " Regula Falsi", &
60 " Pegasus", &
61 " King", &
62 " Anderson Bjoerk", &
63 " Ridder", &
64 " Brent Dekker", &
65 " Newton" /)
66 REAL, PARAMETER, DIMENSION(2,NUM_TESTS) :: bounds = reshape((/ 0.0, 1.2, &
67 0.4, 1.6, &
68 -0.5, 1.9, &
69 -0.5, 0.7, &
70 -1.4, 1.0, &
71 -0.8, 1.6, &
72 -0.5, 1.9, &
73 0.001, 1.201, &
74 -0.9, 1.5, &
75 0.4, 1.0, &
76 -1.2, 0.0, &
77 1.0, 3.4, &
78 0.0, 5.0 /), &
79 (/2,num_tests/))
80 DOUBLE PRECISION, PARAMETER, DIMENSION(NUM_TESTS) :: ref_roots = &
81 (/ 3.9942229171096819451d-01, 8.0413309750366432374d-01, 9.0340766319186021294d-01, &
82 7.7014241346192677110d-02, 2.5920449372984746773d-01, 5.3674166257799978186d-01, &
83 4.4754176206055907112d-01, 1.1111111111111111111d-01, 5.0000003403025908310d-01, &
84 6.7980892150470050192d-01, -3.5938136638046273022d-01, 1.6487212707001281468d-00, &
85 1.0000000000000000000d-00 &
86 /)
87 !--------------------------------------------------------------------------!
88 REAL :: root, xm, dx_rel, dx_acc, plist(1)
89 INTEGER :: i,k, iter, error
90 CHARACTER(LEN=64) :: tap_message
91 LOGICAL :: verbose_results = .false.
92 !--------------------------------------------------------------------------!
93 ! some information
94 WRITE (*,"(A)") "====================================================================="
95 WRITE (*,"(A)") " << Testing root finding algorithms >> "
96 WRITE (*,"(A)") "====================================================================="
97
98
99 tap_plan((num_methods-1)*num_tests)
100
101 DO k=1,num_tests
102 ! plist is an array of real parameters that can be passed to the function func
103 ! so far just a dummy, but it can be used if necessary;
104 plist(1) = 1.0
105 ! reset error code; if GetRoot returns error .ne. 0 something bad happend
106 error = 0
107 ! initial guess: arithmetic mean of boundary values
108 xm = 0.5*(bounds(1,k)+bounds(2,k))
109 print '(A,I2,A,2(ES10.2))', "Test #",k," Search Interval: ",bounds(1,k),bounds(2,k)
110 print '(A)', "---------------------------------------------------------------------"
111 IF (verbose_results) print '(A20,ES27.19)', "reference result", ref_roots(k)
112 DO i=1,num_methods
113 SELECT CASE(i)
114 CASE(1)
115 CALL getroot_bisection(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
116 CASE(2)
117 ! don't use Regula Falsi because of slow convergence in some cases
118 CALL getroot_regulafalsi(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
119 CASE(3)
120 CALL getroot_pegasus(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
121 CASE(4)
122 CALL getroot_king(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
123 CASE(5)
124 CALL getroot_andersonbjoerk(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
125 CASE(6)
126 CALL getroot_ridder(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
127 CASE(7)
128 CALL getroot_brentdekker(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
129 CASE(8)
130 CALL getroot_newton(funcd,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
131 END SELECT
132
133 ! compute the relative error using reference root
134 dx_rel = abs(1.0-root/ref_roots(k))
135 IF (verbose_results) CALL printresults(method_name(i),root,ref_roots(k),error,iter)
136 ! exclude Regular Falsi from check
137 IF (i.NE.2) THEN
138
139 dx_acc = default_accuracy/abs(ref_roots(k))
140 IF (k.EQ.12) THEN
141 ! this test never yields the required accuracy, see \cite engeln2011
142 dx_acc = 1.0e-5
143 END IF
144 WRITE (tap_message,'(A16,A11,ES9.2,A3,ES9.2,A7,I4)') method_name(i), ": dx_rel = ", &
145 dx_rel, " < ", dx_acc, " iter = ", iter
146 tap_check(dx_rel.LE.dx_acc,tap_message)
147
148 END IF
149 END DO
150 print '(A)', "---------------------------------------------------------------------"
151 END DO
152
153 tap_done
154
155CONTAINS
156
157
158 SUBROUTINE printresults(method,root,ref_root,error,iter)
159 IMPLICIT NONE
160 CHARACTER(LEN=*), INTENT(IN) :: method
161 REAL, INTENT(IN) :: root
162 DOUBLE PRECISION, INTENT(IN) :: ref_root
163 INTEGER, INTENT(IN) :: error,iter
164 IF (error.NE.0) THEN
165 print *,trim(geterrormessage(error))
166 ELSE
167 print '(A20,ES23.15,ES9.1,I4)',trim(method),root,abs(1.0-root/ref_root),iter
168 END IF
169 END SUBROUTINE printresults
170
171
172! TAP_PLAN(7)
173!
174!
175!
176! ! Check if the random numbers are in (0,1)
177! ! and if the average is near 0.5
178! TAP_CHECK_CLOSE(r/imax,0.5,1.E-4,"Average close to 0.5.")
179! TAP_CHECK_GE(rmin,0.,"All are bigger (or equal) than 0.")
180! TAP_CHECK_LE(root,1.,"All are smaller (or equal) than 1.")
181! TAP_CHECK_CLOSE(rmin,0.,1.E-4,"Lower limit is close to 0.")
182! TAP_CHECK_CLOSE(rmax,1.,1.E-4,"Upper limit is close to 1.")
183!
184! TAP_CHECK(x.EQ.x0,"SuperKiss64")
185!
186
187END PROGRAM rootstest
188
189
190! function definitions
191PURE SUBROUTINE func(x,fx,plist)
192IMPLICIT NONE
193 COMMON /test_num/ k
194 REAL, INTENT(IN) :: x
195 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
196 REAL, INTENT(OUT) :: fx
197 INTEGER :: k
198 SELECT CASE(k)
199 CASE(1)
200 fx = x*x*(x*x/3. + sqrt(2.)*sin(x)) - sqrt(3.)/18.
201 CASE(2)
202 fx = 11*x**11 - 1.0
203 CASE(3)
204 fx = 35*x**35 - 1.0
205 CASE(4)
206 fx = 2*(x*exp(-9.)-exp(-9*x)) + 1.
207 CASE(5)
208 fx = x*x - (1.-x)**9
209 CASE(6)
210 fx = (x-1.)*exp(-9*x) + x**9
211 CASE(7)
212 fx = x*x + sin(x/9.) - 0.25
213 CASE(8)
214 fx = 0.125*(9.-1./x)
215 CASE(9)
216 fx = tan(x) - x - 0.0463025
217 CASE(10)
218 fx = x*(x+sin(sqrt(75.)*x)) - 0.2
219 CASE(11)
220 fx = x**9 + 1e-4
221 CASE(12)
222 fx = log(x) + 0.5*x*x/exp(1.) -2*x/sqrt(exp(1.)) + 1.
223 CASE DEFAULT
224 fx = x - 1.0
225 END SELECT
226END SUBROUTINE func
227
228
229PURE SUBROUTINE funcd(x,fx,dfx,plist)
230 IMPLICIT NONE
231 COMMON /test_num/ k
232 REAL, INTENT(IN) :: x
233 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
234 REAL, INTENT(OUT) :: fx,dfx
235 INTEGER :: k
236 SELECT CASE(k)
237 CASE(1)
238 fx = x*x*(x*x/3. + sqrt(2.)*sin(x)) - sqrt(3.)/18.
239 dfx= 2*x*(sqrt(2.)*sin(x)+x*x/3.)+x*x*(sqrt(2.)*cos(x)+(2*x)/3.)
240 CASE(2)
241 fx = 11*x**11 - 1.0
242 dfx = 121*x**10
243 CASE(3)
244 fx = 35*x**35 - 1.0
245 dfx= 35*35*x**34 - 1.0
246 CASE(4)
247 fx = 2*(x*exp(-9.)-exp(-9*x)) + 1.
248 dfx= 2*(9*exp(-9*x)+exp(-9.))
249 CASE(5)
250 fx = x*x - (1.-x)**9
251 dfx= 2*x + 9*(1.-x)**8
252 CASE(6)
253 fx = (x-1.)*exp(-9*x) + x**9
254 dfx= (10.-9*x)*exp(-9*x) + 9*x**8
255 CASE(7)
256 fx = x*x + sin(x/9.) - 0.25
257 dfx= 2*x + cos(x/9.)/9.
258 CASE(8)
259 fx = 0.125*(9.-1./x)
260 dfx= 0.125/x**2
261 CASE(9)
262 fx = tan(x) - x - 0.0463025
263 dfx= 1./cos(x)**2 - 1.
264 CASE(10)
265 fx = x*(x+sin(sqrt(75.)*x)) - 0.2
266 dfx= sin(sqrt(75.)*x) + x*(2. + sqrt(75.)*cos(sqrt(75.)*x))
267 CASE(11)
268 fx = x**9 + 1e-4
269 dfx= 9*x**8
270 CASE(12)
271 fx = log(x) + 0.5*x*x/exp(1.) - 2*x/sqrt(exp(1.)) + 1.
272 dfx= 1./x + x/exp(1.) - 2./sqrt(exp(1.))
273 CASE DEFAULT
274 fx = x - 1.0
275 dfx= 1.
276 END SELECT
277END SUBROUTINE funcd
278
pure subroutine funcd(y, fy, dfy, plist)
Definition: bondi2d.f90:416
root finding subroutines
Definition: roots.f90:45
subroutine, public getroot_brentdekker(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:485
subroutine, public getroot_bisection(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:519
subroutine, public getroot_regulafalsi(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:315
subroutine, public getroot_king(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:383
subroutine, public getroot_ridder(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:451
subroutine, public getroot_andersonbjoerk(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:417
subroutine, public getroot_newton(funcd, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:247
subroutine, public getroot_pegasus(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:349
character(len=64) function, public geterrormessage(error)
Definition: roots.f90:160
real, parameter, public default_accuracy
Definition: roots.f90:67
pure subroutine func(x, fx, plist)
Definition: rootstest.f90:192
program rootstest
Definition: rootstest.f90:32
subroutine printresults(method, root, ref_root, error, iter)
Definition: rootstest.f90:159