Function to be analyzed. f(ax) and f(bx) should have opposite signs, or lowerBound and/or upperBound should be defined to perform initial interval extension.
Defines an early termination condition. Receives the current upper and lower bounds on the root. The delegate must return true when these bounds are acceptable. If this function always returns false or it is null, full machine precision will be achieved.
Left inner bound of initial range of f known to contain the root.
Right inner bound of initial range of f known to contain the root. Can be equal to ax.
Value of f(ax) (optional).
Value of f(bx) (optional).
Lower outer bound for interval extension (optional)
Upper outer bound for interval extension (optional)
Appr. maximum allowed number of function calls (inluciding function calls for grid).
Number of nodes in the left side and right side regular grids (optional). If it equals to 0 (default), then the algorithm uses ieeeMean for searching. The algorithm performs searching if sgn(fax) equals to sgn(fbx) and at least one outer bound allows to extend the inner initial range.
import mir.math.common: log, exp, fabs; auto logRoot = findRoot!log(0, double.infinity).validate.x; assert(logRoot == 1); auto shift = 1; auto expm1Root = findRoot!(x => exp(x) - shift) (-double.infinity, double.infinity).validate.x; assert(expm1Root == 0); auto approxLogRoot = findRoot!(log, (a, b) => fabs(a - b) < 1e-5)(0, double.infinity).validate.x; assert(fabs(approxLogRoot - 1) < 1e-5);
With adaptive bounds
import mir.math.common: log, exp, fabs; auto logRoot = findRoot!log( 10, 10, // assume we have one initial point double.nan, double.nan, // fa, fb aren't provided by user -double.infinity, double.infinity, // all space is available for the bounds extension. ).validate.x; assert(logRoot == 1); auto shift = 1; alias expm1Fun = (double x) => exp(x) - shift; auto expm1RootRet = findRoot!expm1Fun ( 11, 10, // reversed order for interval is always OK expm1Fun(11), expm1Fun(10), // the order must be the same as bounds 0, double.infinity, // space for the bounds extension. ).validate; assert(expm1Fun(expm1RootRet.x) == 0); auto approxLogRoot = findRoot!(log, (a, b) => fabs(a - b) < 1e-5)( -1e10, +1e10, double.nan, double.nan, 0, double.infinity, ).validate.x; assert(fabs(approxLogRoot - 1) < 1e-5);
ditto
1 import core.stdc.tgmath: atan; 2 import mir.math; 3 import std.meta: AliasSeq; 4 5 const double[2][3] boundaries = [ 6 [0.4, 0.6], 7 [1.4, 1.6], 8 [0.1, 2.1]]; 9 10 enum root = 1.0; 11 12 foreach(fun; AliasSeq!( 13 (double x) => x ^^ 2 - root, 14 (double x) => root - x ^^ 2, 15 (double x) => atan(x - root), 16 )) 17 { 18 foreach(ref bounds; boundaries) 19 { 20 auto result = findRoot!fun( 21 bounds[0], bounds[1], 22 double.nan, double.nan, // f(a) and f(b) not provided 23 -double.max, double.max, // user provided outer bounds 24 ); 25 assert(result.validate.x == root); 26 } 27 } 28 29 foreach(ref bounds; boundaries) 30 { 31 auto result = findRoot!(x => sin(x - root))( 32 bounds[0], bounds[1], 33 double.nan, double.nan, // f(a) and f(b) not provided 34 -10, 10, // user provided outer bounds 35 100, // max iterations, 36 10, // steps for grid 37 ); 38 assert(result.validate.x == root); 39 } 40 // single initial point, grid, positive direction 41 auto result = findRoot!((double x ) => sin(x - root))( 42 0.1, 0.1, // initial point, a = b 43 double.nan, double.nan, // f(a) = f(b) not provided 44 -100.0, 100.0, // user provided outer bounds 45 150, // max iterations, 46 100, // number of steps for grid 47 ); 48 assert(result.validate.x == root); 49 50 // single initial point, grid, negative direction 51 result = findRoot!((double x ) => sin(x - root))( 52 0.1, 0.1, // initial point a = b 53 double.nan, double.nan, // f(a) = f(b) not provided 54 100.0, -100.0, // user provided outer bounds, reversed order 55 150, // max iterations, 56 100, // number of steps for grid 57 ); 58 assert(result.validate.x == double(root - PI)); // Left side root!
With adaptive bounds and single initial point. Reverse outer bound order controls first step direction in case of f(a) == f(b).
enum root = 1.0; // roots are +/- `root` alias fun = (double x) => x * x - root; double lowerBound = -10.0; double upperBound = 10.0; assert( findRoot!fun( 0, 0, // initial interval double.nan, double.nan, lowerBound, upperBound, // positive direction has higher priority ).validate.x == root ); assert( findRoot!fun( 0, 0, // initial interval double.nan, double.nan, upperBound, lowerBound, // reversed order ).validate.x == -root // other root );
Find root of a real function f(x) by bracketing, allowing the termination condition to be specified.
Given a function f and a range [a .. b] such that f(a) and f(b) have opposite signs or at least one of them equals ±0, returns the value of x in the range which is closest to a root of f(x). If f(x) has more than one root in the range, one will be chosen arbitrarily. If f(x) returns NaN, NaN will be returned; otherwise, this algorithm is guaranteed to succeed.
Uses an algorithm based on TOMS748, which uses inverse cubic interpolation whenever possible, otherwise reverting to parabolic or secant interpolation. Compared to TOMS748, this implementation improves worst-case performance by a factor of more than 100, and typical performance by a factor of 2. For 80-bit reals, most problems require 8 to 15 calls to f(x) to achieve full machine precision. The worst-case performance (pathological cases) is approximately twice the number of bits.
References: "On Enclosing Simple Roots of Nonlinear Equations", G. Alefeld, F.A. Potra, Yixun Shi, Mathematics of Computation 61, pp733-744 (1993). Fortran code available from www.netlib.org as algorithm TOMS478.