ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
pm_mathRoot::setRoot Interface Reference

Return a root of a specified continuous real-valued one-dimensional mathematical function such that \(f(\mathrm{root}) = 0\) with the user-specified or the default root-finding method.
More...

Detailed Description

Return a root of a specified continuous real-valued one-dimensional mathematical function such that \(f(\mathrm{root}) = 0\) with the user-specified or the default root-finding method.

See the documentation of pm_mathRoot for details of the various root-finding methods.

Note
Which root-finding method should I use among all available?
The brent method is the recommended algorithm of choice for general one-dimensional root-finding problems without knowledge of function derivative (gradient).
The Newton-Raphson method is the recommended algorithm of choice for general one-dimensional root-finding problems with knowledge of function derivative (gradient).
Parameters
[in]method: The input scalar constant that can be one of the following:
  1. The constant false or an object of type false_type, signifying the use of the False-Position root-finding method within the algorithm.
  2. The constant bisection or an object of type bisection_type, signifying the use of the Bisection root-finding method within the algorithm.
  3. The constant secant or an object of type secant_type, signifying the use of the Secant root-finding method within the algorithm.
  4. The constant brent or an object of type brent_type, signifying the use of the Brent root-finding method within the algorithm.
  5. The constant ridders or an object of type ridders_type, signifying the use of the Ridders root-finding method within the algorithm.
  6. The constant toms748 or an object of type toms748_type, signifying the use of the TOMS748 root-finding method within the algorithm.
  7. The constant newton or an object of type newton_type, signifying the use of the Newton root-finding method within the algorithm.
  8. The constant halley or an object of type halley_type, signifying the use of the Halley root-finding method within the algorithm.
  9. The constant schroder or an object of type schroder_type, signifying the use of the Schroder root-finding method within the algorithm.
[in]getFunc: The external user-specified function whose interface depends on the specified value for the input argument method.
  1. If the specified method is any of the following,
    1. the constant brent or an object of type brent_type,
    2. the constant false or an object of type false_type,
    3. the constant secant or an object of type secant_type,
    4. the constant ridders or an object of type ridders_type,
    5. the constant toms748 or an object of type toms748_type,
    6. the constant bisection or an object of type bisection_type,
    none of which require the derivative of the target function to find the function root, then getFunc() must take a single scalar input argument x of the same type and kind as the output argument root (below).
    On output, getFunc() must return a scalar of the same type and kind as the function input argument x, containing the value of the target function evaluated at the specified input x.
    The following illustrates the generic interface of getFunc() for the above values of method,
    function getFunc(x) result(func)
    real(RKG) , intent(in) :: x
    real(RKG) :: func
    end function
    where RKG refers to the kind of the output argument root.
  2. If the specified method is any of the following,
    1. the constant newton or an object of type newton_type,
    2. the constant halley or an object of type halley_type,
    3. the constant schroder or an object of type schroder_type,
    all of which either require the first (Newton) or higher derivatives (Halley/Schroder) of the target function to find the function root, then getFunc() must take two arguments:
    1. A scalar input argument x of the same type and kind as the output argument root (below).
    2. A scalar input argument order of type integer of default kind IK.
    On output, getFunc() must return a scalar of the same type and kind as the function input argument x, containing the derivative of the specified input order of the target function evaluated at the specified input x:
    1. An input order value of 0 must yield the target function value at the input x.
    2. An input order value of 1 must yield the target function first derivative at the input x.
    3. An input order value of 2 must yield the target function second derivative at the input x.
    The following illustrates the generic interface of getFunc() for the above values of method,
    function getFunc(x, order) result(func)
    real(RKG) , intent(in) :: x
    integer(IK) , intent(in) :: order
    real(RKG) :: func
    end function
    where RKG refers to the kind of the output argument root.
[in,out]root: The input/output scalar of type real of kind any supported by the processor (e.g., RK, RK32, RK64, or RK128).
On input,
  1. If the specified method is of type newton_type or halley_type, and if lb < root .and. root <ub, then the specified input value for root will be used as the best guess initial value for the function root within the bracket specified by the input arguments lb and ub.
    Otherwise the center of the specifie dinput bracket [lb, ub] wil be used as the initial best-guess.
  2. Any input value for root is ignored for all other root-finding methods.
On output, root contains the best approximation to the root of the user-specified function getFunc() within the specified search interval [lb, ub] and tolerance threshold abstol.
[in]lb: The input scalar of same type and kind as the output root, representing the lower bound of the bracket (search) interval.
[in]ub: The input scalar of same type and kind as the output root, representing the upper bound of the bracket (search) interval.
[in]abstol: The input scalar of same type and kind as the output root, representing the absolute tolerance used as the stopping criterion of the search.
The iterations of the specified input method continue until the search interval becomes smaller than abstol in absolute units.
Care must be taken for specifying a reasonable value for abstol (see the warnings below).
If no suitable value for abstol is known a priori, try abstol = epsilon(0._RKG)**.8 * (abs(lb) + abs(ub)) where RKG refers to the kind of the output argument root.
[out]neval: The output scalar argument of type integer of default kind IK, containing the total number of getFunc() function calls made by the algorithm.
  1. A positive output neval implies the successful convergence of the algorithm to the function root after +neval function evaluations.
  2. A negative output neval implies the failure of convergence of the algorithm to the function root after -neval function evaluations.
  3. A zero output neval only occurs if the root lies at the search bracket boundaries specified by the input argument lb or ub.
Convergence failures rarely occur. If they do, setting the input arguments abstol and/or niter to larger values may resolve the failure.
[in]niter: The input scalar of type integer of default kind IK, representing the maximum number of iterations allowed for the specified method to find the root.
The default number of steps niter within the algorithm is a compile-time constant that depends on the kind of the real input arguments.
(optional, default = ceiling(2 * precision(lb) / log10(2.)))


Possible calling interfaces

use pm_mathRoot, only: setRoot
call setRoot(method, getFunc, root, lb, ub, lf, uf, abstol, neval)
call setRoot(method, getFunc, root, lb, ub, lf, uf, abstol, neval, niter)
Return a root of a specified continuous real-valued one-dimensional mathematical function such that ...
This module contains classes and procedures for computing the roots of one-dimensional continuous mat...
Warning
The condition lb < ub must hold for the corresponding input arguments.
The condition 0._RKG < abstol must hold for the corresponding input arguments.
The condition abstol < (ub - lb) must hold for the corresponding input arguments.
The condition abs(getFunc(lb), lf) < abstol must hold for the corresponding input arguments.
The condition abs(getFunc(ub), uf) < abstol must hold for the corresponding input arguments.
The condition sign(lf, 1.) /= sign(uf, 1.) must hold for the corresponding input arguments.
The condition 0 < niter must hold for the corresponding input arguments.
These conditions are verified only if the library is built with the preprocessor macro CHECK_ENABLED=1.
It is crucial to keep in mind that computers use a fixed number of binary digits to represent floating-point numbers.
While the user-specified function might analytically pass through zero, it is possible that its computed value is never zero, for any floating-point argument.
One must also decide on what accuracy on the root is attainable.
For example, convergence to within \(10^{−6}\) in absolute value is reasonable when the root lies near \(1\), but unachievable if the root lies near an extremely large number near floating-point overflow.
Remarks
The procedures under discussion are impure.
The procedures under discussion are recursive.
See also
getRoot
setRoot
getPolyRoot
setPolyRoot


Example usage

1program example
2
3 use pm_kind, only: SK, IK, LK
4 use pm_kind, only: RKG => RKH ! all processor kinds are supported.
6 use pm_io, only: display_type
7
8 implicit none
9
10 real(RKG) :: root
11 real(RKG) :: lb, ub
12 real(RKG) :: lf, uf
13 logical(LK) :: failed
14 integer(IK) :: neval
15 type(display_type) :: disp
16 disp = display_type(file = "main.out.F90")
17
18 call disp%skip()
19 call disp%show("lb = 2._RKG; lf = getSin(lb)")
20 lb = 2._RKG; lf = getSin(lb)
21 call disp%show("ub = 4._RKG; uf = getSin(ub)")
22 ub = 4._RKG; uf = getSin(ub)
23 call disp%show("[lf, uf]")
24 call disp%show( [lf, uf] )
25 call disp%show("call setRoot(false, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
26 call setRoot(false, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
27 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
28 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
29 call disp%show("[root, getSin(root)]")
30 call disp%show( [root, getSin(root)] )
31 call disp%show("neval")
32 call disp%show( neval )
33 call disp%skip()
34 call disp%show("call setRoot(bisection, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
35 call setRoot(bisection, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
36 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
37 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
38 call disp%show("[root, getSin(root)]")
39 call disp%show( [root, getSin(root)] )
40 call disp%show("neval")
41 call disp%show( neval )
42 call disp%skip()
43 call disp%show("call setRoot(secant, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
44 call setRoot(secant, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
45 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
46 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
47 call disp%show("[root, getSin(root)]")
48 call disp%show( [root, getSin(root)] )
49 call disp%show("neval")
50 call disp%show( neval )
51 call disp%skip()
52 call disp%show("call setRoot(brent, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
53 call setRoot(brent, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
54 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
55 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
56 call disp%show("[root, getSin(root)]")
57 call disp%show( [root, getSin(root)] )
58 call disp%show("neval")
59 call disp%show( neval )
60 call disp%skip()
61 call disp%show("call setRoot(ridders, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
62 call setRoot(ridders, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
63 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
64 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
65 call disp%show("[root, getSin(root)]")
66 call disp%show( [root, getSin(root)] )
67 call disp%show("neval")
68 call disp%show( neval )
69 call disp%skip()
70 call disp%show("call setRoot(toms748, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
71 call setRoot(toms748, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
72 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
73 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
74 call disp%show("[root, getSin(root)]")
75 call disp%show( [root, getSin(root)] )
76 call disp%show("neval")
77 call disp%show( neval )
78 call disp%skip()
79 call disp%show("call setRoot(newton, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
80 call setRoot(newton, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
81 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
82 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
83 call disp%show("[root, getSin(root)]")
84 call disp%show( [root, getSin(root)] )
85 call disp%show("neval")
86 call disp%show( neval )
87 call disp%skip()
88 call disp%show("call setRoot(halley, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
89 call setRoot(halley, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
90 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
91 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
92 call disp%show("[root, getSin(root)]")
93 call disp%show( [root, getSin(root)] )
94 call disp%show("neval")
95 call disp%show( neval )
96 call disp%skip()
97 call disp%show("call setRoot(schroder, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
98 call setRoot(schroder, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
99 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
100 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
101 call disp%show("[root, getSin(root)]")
102 call disp%show( [root, getSin(root)] )
103 call disp%show("neval")
104 call disp%show( neval )
105 call disp%skip()
106
107
108
109 call disp%skip()
110 call disp%show("lb = 1._RKG; lf = getCos(lb)")
111 lb = 1._RKG; lf = getCos(lb)
112 call disp%show("ub = 3._RKG; uf = getCos(ub)")
113 ub = 3._RKG; uf = getCos(ub)
114 call disp%show("[lf, uf]")
115 call disp%show( [lf, uf] )
116 call disp%show("call setRoot(false, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
117 call setRoot(false, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
118 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
119 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
120 call disp%show("[root, getCos(root)]")
121 call disp%show( [root, getCos(root)] )
122 call disp%show("neval")
123 call disp%show( neval )
124 call disp%skip()
125 call disp%show("call setRoot(bisection, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
126 call setRoot(bisection, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
127 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
128 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
129 call disp%show("[root, getCos(root)]")
130 call disp%show( [root, getCos(root)] )
131 call disp%show("neval")
132 call disp%show( neval )
133 call disp%skip()
134 call disp%show("call setRoot(secant, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
135 call setRoot(secant, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
136 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
137 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
138 call disp%show("[root, getCos(root)]")
139 call disp%show( [root, getCos(root)] )
140 call disp%show("neval")
141 call disp%show( neval )
142 call disp%skip()
143 call disp%show("call setRoot(brent, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
144 call setRoot(brent, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
145 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
146 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
147 call disp%show("[root, getCos(root)]")
148 call disp%show( [root, getCos(root)] )
149 call disp%show("neval")
150 call disp%show( neval )
151 call disp%skip()
152 call disp%show("call setRoot(ridders, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
153 call setRoot(ridders, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
154 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
155 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
156 call disp%show("[root, getCos(root)]")
157 call disp%show( [root, getCos(root)] )
158 call disp%show("neval")
159 call disp%show( neval )
160 call disp%skip()
161 call disp%show("call setRoot(toms748, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
162 call setRoot(toms748, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
163 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
164 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
165 call disp%show("[root, getCos(root)]")
166 call disp%show( [root, getCos(root)] )
167 call disp%show("neval")
168 call disp%show( neval )
169 call disp%skip()
170 call disp%show("call setRoot(newton, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
171 call setRoot(newton, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
172 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
173 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
174 call disp%show("[root, getCos(root)]")
175 call disp%show( [root, getCos(root)] )
176 call disp%show("neval")
177 call disp%show( neval )
178 call disp%skip()
179 call disp%show("call setRoot(halley, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
180 call setRoot(halley, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
181 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
182 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
183 call disp%show("[root, getCos(root)]")
184 call disp%show( [root, getCos(root)] )
185 call disp%show("neval")
186 call disp%show( neval )
187 call disp%skip()
188 call disp%show("call setRoot(schroder, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
189 call setRoot(schroder, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
190 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
191 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
192 call disp%show("[root, getCos(root)]")
193 call disp%show( [root, getCos(root)] )
194 call disp%show("neval")
195 call disp%show( neval )
196 call disp%skip()
197
198
199
200 call disp%skip()
201 call disp%show("lb = -1._RKG; lf = getQuad(lb)")
202 lb = -1._RKG; lf = getQuad(lb)
203 call disp%show("ub = +4._RKG; uf = getQuad(ub)")
204 ub = +4._RKG; uf = getQuad(ub)
205 call disp%show("[lf, uf]")
206 call disp%show( [lf, uf] )
207 call disp%show("call setRoot(false, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
208 call setRoot(false, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
209 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
210 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
211 call disp%show("[root, getQuad(root)]")
212 call disp%show( [root, getQuad(root)] )
213 call disp%show("neval")
214 call disp%show( neval )
215 call disp%skip()
216 call disp%show("call setRoot(bisection, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
217 call setRoot(bisection, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
218 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
219 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
220 call disp%show("[root, getQuad(root)]")
221 call disp%show( [root, getQuad(root)] )
222 call disp%show("neval")
223 call disp%show( neval )
224 call disp%skip()
225 call disp%show("call setRoot(secant, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
226 call setRoot(secant, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
227 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
228 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
229 call disp%show("[root, getQuad(root)]")
230 call disp%show( [root, getQuad(root)] )
231 call disp%show("neval")
232 call disp%show( neval )
233 call disp%skip()
234 call disp%show("call setRoot(brent, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
235 call setRoot(brent, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
236 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
237 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
238 call disp%show("[root, getQuad(root)]")
239 call disp%show( [root, getQuad(root)] )
240 call disp%show("neval")
241 call disp%show( neval )
242 call disp%skip()
243 call disp%show("call setRoot(ridders, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
244 call setRoot(ridders, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
245 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
246 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
247 call disp%show("[root, getQuad(root)]")
248 call disp%show( [root, getQuad(root)] )
249 call disp%show("neval")
250 call disp%show( neval )
251 call disp%skip()
252 call disp%show("call setRoot(toms748, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
253 call setRoot(toms748, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
254 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
255 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
256 call disp%show("[root, getQuad(root)]")
257 call disp%show( [root, getQuad(root)] )
258 call disp%show("neval")
259 call disp%show( neval )
260 call disp%skip()
261 call disp%show("call setRoot(newton, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
262 call setRoot(newton, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
263 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
264 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
265 call disp%show("[root, getQuad(root)]")
266 call disp%show( [root, getQuad(root)] )
267 call disp%show("neval")
268 call disp%show( neval )
269 call disp%skip()
270 call disp%show("call setRoot(halley, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
271 call setRoot(halley, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
272 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
273 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
274 call disp%show("[root, getQuad(root)]")
275 call disp%show( [root, getQuad(root)] )
276 call disp%show("neval")
277 call disp%show( neval )
278 call disp%skip()
279 call disp%show("call setRoot(schroder, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)")
280 call setRoot(schroder, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
281 call disp%show("if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'")
282 if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
283 call disp%show("[root, getQuad(root)]")
284 call disp%show( [root, getQuad(root)] )
285 call disp%show("neval")
286 call disp%show( neval )
287 call disp%skip()
288
289contains
290
291 pure function getSin(x) result(func)
292 real(RKG), intent(in) :: x
293 real(RKG) :: func
294 func = sin(x)
295 end function
296
297 pure function getSinDiff(x, order) result(func)
298 integer(IK), intent(in) :: order
299 real(RKG), intent(in) :: x
300 real(RKG) :: func
301 if (order == 0) func = +getSin(x)
302 if (order == 1) func = +cos(x)
303 if (order == 2) func = -sin(x)
304 end function
305
306 pure function getCos(x) result(func)
307 real(RKG), intent(in) :: x
308 real(RKG) :: func
309 func = cos(x)
310 end function
311
312 pure function getCosDiff(x, order) result(func)
313 integer(IK), intent(in) :: order
314 real(RKG), intent(in) :: x
315 real(RKG) :: func
316 if (order == 0) func = +getCos(x)
317 if (order == 1) func = -sin(x)
318 if (order == 2) func = -cos(x)
319 end function
320
321 pure function getQuad(x) result(func)
322 real(RKG), intent(in) :: x
323 real(RKG) :: func
324 func = x * (x - 1._RKG) * (x - 2._RKG)
325 end function
326
327 pure function getQuadDiff(x, order) result(func)
328 integer(IK), intent(in) :: order
329 real(RKG), intent(in) :: x
330 real(RKG) :: func
331 if (order == 0) func = getQuad(x)
332 if (order == 1) func = 3._RKG * x**2 - 6._RKG * x + 2._RKG
333 if (order == 2) func = 6._RKG * x - 6._RKG
334 end function
335
336end program example
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
integer, parameter RKH
The scalar integer constant of intrinsic default kind, representing the highest-precision real kind t...
Definition: pm_kind.F90:858
type(brent_type), parameter brent
This is a scalar parameter object of type brent_type that is exclusively used to signify the use of B...
type(ridders_type), parameter ridders
This is a scalar parameter object of type ridders_type that is exclusively used to signify the use of...
type(toms748_type), parameter toms748
This is a scalar parameter object of type toms748_type that is exclusively used to signify the use of...
type(newton_type), parameter newton
This is a scalar parameter object of type newton_type that is exclusively used to signify the use of ...
type(bisection_type), parameter bisection
This is a scalar parameter object of type bisection_type that is exclusively used to signify the use ...
type(false_type), parameter false
This is a scalar parameter object of type false_type that is exclusively used to signify the use of F...
type(halley_type), parameter halley
This is a scalar parameter object of type halley_type that is exclusively used to signify the use of ...
type(schroder_type), parameter schroder
This is a scalar parameter object of type schroder_type that is exclusively used to signify the use o...
type(secant_type), parameter secant
This is a scalar parameter object of type secant_type that is exclusively used to signify the use of ...
Generate and return an object of type display_type.
Definition: pm_io.F90:10282

Example Unix compile command via Intel ifort compiler
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example output
1
2lb = 2._RKG; lf = getSin(lb)
3ub = 4._RKG; uf = getSin(ub)
4[lf, uf]
5+0.909297426825681695396019865911744868, -0.756802495307928251372639094511829109
6call setRoot(false, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
7if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
8[root, getSin(root)]
9+3.14159265358979323846264338327950280, +0.867181013012378102479704402604335225E-34
10neval
11+7
12
13call setRoot(bisection, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
14if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
15[root, getSin(root)]
16+3.14159265358979323846264376397729529, -0.380697792403865321819264032602920990E-24
17neval
18+82
19
20call setRoot(secant, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
21if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
22[root, getSin(root)]
23+3.14159265358979323846264338327950280, +0.867181013012378102479704402604335225E-34
24neval
25+7
26
27call setRoot(brent, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
28if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
29[root, getSin(root)]
30+3.14159265358979323846264338327950280, +0.867181013012378102479704402604335225E-34
31neval
32+7
33
34call setRoot(ridders, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
35if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
36[root, getSin(root)]
37+3.14159265358979323846264338327950280, +0.867181013012378102479704402604335225E-34
38neval
39+20
40
41call setRoot(toms748, getSin, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
42if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
43[root, getSin(root)]
44+3.14159265358979323846264338327950280, +0.867181013012378102479704402604335225E-34
45neval
46+7
47
48call setRoot(newton, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
49if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
50[root, getSin(root)]
51+3.14159265358979323846264338327950280, +0.867181013012378102479704402604335225E-34
52neval
53+10
54
55call setRoot(halley, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
56if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
57[root, getSin(root)]
58+3.14159265358979323846264338327950280, +0.867181013012378102479704402604335225E-34
59neval
60+3
61
62call setRoot(schroder, getSinDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
63if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
64[root, getSin(root)]
65+3.14159265358979323846264338327950280, +0.867181013012378102479704402604335225E-34
66neval
67+3
68
69
70lb = 1._RKG; lf = getCos(lb)
71ub = 3._RKG; uf = getCos(ub)
72[lf, uf]
73+0.540302305868139717400936607442976558, -0.989992496600445457271572794731261336
74call setRoot(false, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
75if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
76[root, getCos(root)]
77+1.57079632679489661923132169163975140, +0.433590506506189051239852201302167613E-34
78neval
79+7
80
81call setRoot(bisection, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
82if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
83[root, getCos(root)]
84+1.57079632679489661923132270916926020, -0.101752950875496033578104070837142335E-23
85neval
86+82
87
88call setRoot(secant, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
89if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
90[root, getCos(root)]
91+1.57079632679489661923132169163975140, +0.433590506506189051239852201302167613E-34
92neval
93+7
94
95call setRoot(brent, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
96if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
97[root, getCos(root)]
98+1.57079632679489661923132169163975140, +0.433590506506189051239852201302167613E-34
99neval
100+7
101
102call setRoot(ridders, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
103if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
104[root, getCos(root)]
105+1.57079632679489661923132169164015931, -0.407868603170565934772132143019357364E-30
106neval
107+20
108
109call setRoot(toms748, getCos, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
110if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
111[root, getCos(root)]
112+1.57079632679489661923132169163975140, +0.433590506506189051239852201302167613E-34
113neval
114+8
115
116call setRoot(newton, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
117if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
118[root, getCos(root)]
119+1.57079632679489661923132169163975140, +0.433590506506189051239852201302167613E-34
120neval
121+10
122
123call setRoot(halley, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
124if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
125[root, getCos(root)]
126+1.57079632679489661923132169163975140, +0.433590506506189051239852201302167613E-34
127neval
128+3
129
130call setRoot(schroder, getCosDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
131if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
132[root, getCos(root)]
133+1.57079632679489661923132169163975140, +0.433590506506189051239852201302167613E-34
134neval
135+3
136
137
138lb = -1._RKG; lf = getQuad(lb)
139ub = +4._RKG; uf = getQuad(ub)
140[lf, uf]
141-6.00000000000000000000000000000000000, +24.0000000000000000000000000000000000
142call setRoot(false, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
143if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
144[root, getQuad(root)]
145+0.00000000000000000000000000000000000, +0.00000000000000000000000000000000000
146neval
147+1
148
149call setRoot(bisection, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
150if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
151[root, getQuad(root)]
152+1.99999999999999999999999958640969372, -0.827180612553027674871408178899138516E-24
153neval
154+83
155
156call setRoot(secant, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
157if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
158[root, getQuad(root)]
159+0.00000000000000000000000000000000000, +0.00000000000000000000000000000000000
160neval
161+1
162
163call setRoot(brent, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
164if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
165[root, getQuad(root)]
166+0.00000000000000000000000000000000000, +0.00000000000000000000000000000000000
167neval
168+1
169
170call setRoot(ridders, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
171if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
172[root, getQuad(root)]
173+2.00000000000000000000000000000000000, +0.00000000000000000000000000000000000
174neval
175+22
176
177call setRoot(toms748, getQuad, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
178if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
179[root, getQuad(root)]
180+2.00000000000000000000000000000000000, +0.00000000000000000000000000000000000
181neval
182+2
183
184call setRoot(newton, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
185if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
186[root, getQuad(root)]
187+2.00000000000000000000000000000000000, +0.00000000000000000000000000000000000
188neval
189+2
190
191call setRoot(halley, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
192if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
193[root, getQuad(root)]
194+2.00000000000000000000000000000000000, +0.00000000000000000000000000000000000
195neval
196+3
197
198call setRoot(schroder, getQuadDiff, root, lb, ub, lf, uf, abstol = epsilon(0._RKG)**.7, neval = neval)
199if (neval < 0_IK) error stop 'The root-finding algorithm failed to converge.'
200[root, getQuad(root)]
201+2.00000000000000000000000000000000000, +0.00000000000000000000000000000000000
202neval
203+3
204
205
Test:
test_pm_mathRoot
Todo:
High Priority: The current implementation of the Schroder method does not converge for certain problems that the Halley method converges readily.
This may be due to the nature of the Schroder update.
Regardless, this should be further investigated.
For now, Schroder is diverted to Halley.


Final Remarks


If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan

Definition at line 2675 of file pm_mathRoot.F90.


The documentation for this interface was generated from the following file: