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: The constant false or an object of type false_type, signifying the use of the False-Position root-finding method within the algorithm. The constant bisection or an object of type bisection_type, signifying the use of the Bisection root-finding method within the algorithm. The constant secant or an object of type secant_type, signifying the use of the Secant root-finding method within the algorithm. The constant brent or an object of type brent_type, signifying the use of the Brent root-finding method within the algorithm. The constant ridders or an object of type ridders_type, signifying the use of the Ridders root-finding method within the algorithm. The constant toms748 or an object of type toms748_type, signifying the use of the TOMS748 root-finding method within the algorithm. The constant newton or an object of type newton_type, signifying the use of the Newton root-finding method within the algorithm. The constant halley or an object of type halley_type, signifying the use of the Halley root-finding method within the algorithm. 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`. If the specified `method` is any of the following, the constant brent or an object of type brent_type, the constant false or an object of type false_type, the constant secant or an object of type secant_type, the constant ridders or an object of type ridders_type, the constant toms748 or an object of type toms748_type, 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`. If the specified `method` is any of the following, the constant newton or an object of type newton_type, the constant halley or an object of type halley_type, 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: A scalar input argument `x` of the same type and kind as the output argument `root` (below). 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`: An input `order` value of `0` must yield the target function value at the input `x`. An input `order` value of `1` must yield the target function first derivative at the input `x`. 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, 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. 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. A positive output `neval` implies the successful convergence of the algorithm to the function root after `+neval` function evaluations. A negative output `neval` implies the failure of convergence of the algorithm to the function root after `-neval` function evaluations. 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)

pm_mathRoot::setRoot
Return a root of a specified continuous real-valued one-dimensional mathematical function such that ...
Definition: pm_mathRoot.F90:2675

pm_mathRoot
This module contains classes and procedures for computing the roots of one-dimensional continuous mat...
Definition: pm_mathRoot.F90:454

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.

5 use pm_mathRoot, only: setRoot, false, bisection, secant, brent, ridders, toms748, newton, halley, schroder

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

pm_io::show
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726

pm_io::skip
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508

pm_io
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252

pm_io::disp
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393

pm_kind
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268

pm_kind::LK
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541

pm_kind::IK
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

pm_kind::SK
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

pm_kind::RKH
integer, parameter RKH
The scalar integer constant of intrinsic default kind, representing the highest-precision real kind t...
Definition: pm_kind.F90:858

pm_mathRoot::brent
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...
Definition: pm_mathRoot.F90:709

pm_mathRoot::ridders
type(ridders_type), parameter ridders
This is a scalar parameter object of type ridders_type that is exclusively used to signify the use of...
Definition: pm_mathRoot.F90:1325

pm_mathRoot::toms748
type(toms748_type), parameter toms748
This is a scalar parameter object of type toms748_type that is exclusively used to signify the use of...
Definition: pm_mathRoot.F90:797

pm_mathRoot::newton
type(newton_type), parameter newton
This is a scalar parameter object of type newton_type that is exclusively used to signify the use of ...
Definition: pm_mathRoot.F90:1061

pm_mathRoot::bisection
type(bisection_type), parameter bisection
This is a scalar parameter object of type bisection_type that is exclusively used to signify the use ...
Definition: pm_mathRoot.F90:1413

pm_mathRoot::false
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...
Definition: pm_mathRoot.F90:885

pm_mathRoot::halley
type(halley_type), parameter halley
This is a scalar parameter object of type halley_type that is exclusively used to signify the use of ...
Definition: pm_mathRoot.F90:1149

pm_mathRoot::schroder
type(schroder_type), parameter schroder
This is a scalar parameter object of type schroder_type that is exclusively used to signify the use o...
Definition: pm_mathRoot.F90:1237

pm_mathRoot::secant
type(secant_type), parameter secant
This is a scalar parameter object of type secant_type that is exclusively used to signify the use of ...
Definition: pm_mathRoot.F90:973

pm_io::display_type
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.

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.
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.

Copyright: Computational Data Science Lab

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:

src/fortran/main/pm_mathRoot.F90