blob: 61ff12c69c509c0f5eb8040ca0b87855d2e7baa7 [file] [log] [blame]
 // Derived from SciPy's special/c_misc/fsolve.c and special/c_misc/misc.h // https://github.com/scipy/scipy/blob/master/scipy/special/c_misc/fsolve.c // https://github.com/scipy/scipy/blob/master/scipy/special/c_misc/misc.h // Copyright ©2017 The gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package mathext import "math" type objectiveFunc func(float64, []float64) float64 type fSolveResult uint8 const ( // An exact solution was found, in which case the first point on the // interval is the value fSolveExact fSolveResult = iota + 1 // Interval width is less than the tolerance fSolveConverged // Root-finding didn't converge in a set number of iterations fSolveMaxIterations ) const ( machEp = 1.0 / (1 << 53) ) // falsePosition uses a combination of bisection and false position to find a // root of a function within a given interval. This is guaranteed to converge, // and always keeps a bounding interval, unlike Newton's method. Inputs are: // x1, x2: initial bounding interval // f1, f2: value of f() at x1 and x2 // absErr, relErr: absolute and relative errors on the bounding interval // bisectTil: if > 0.0, perform bisection until the width of the bounding // interval is less than this // f, fExtra: function to find root of is f(x, fExtra) // Returns: // result: whether an exact root was found, the process converged to a // bounding interval small than the required error, or the max number // of iterations was hit // bestX: best root approximation // bestF: function value at bestX // errEst: error estimation func falsePosition(x1, x2, f1, f2, absErr, relErr, bisectTil float64, f objectiveFunc, fExtra []float64) (fSolveResult, float64, float64, float64) { // The false position steps are either unmodified, or modified with the // Anderson-Bjorck method as appropriate. Theoretically, this has a "speed of // convergence" of 1.7 (bisection is 1, Newton is 2). // Note that this routine was designed initially to work with gammaincinv, so // it may not be tuned right for other problems. Don't use it blindly. if f1*f2 >= 0 { panic("Initial interval is not a bounding interval") } const ( maxIterations = 100 bisectIter = 4 bisectWidth = 4.0 ) const ( bisect = iota + 1 falseP ) var state uint8 if bisectTil > 0 { state = bisect } else { state = falseP } gamma := 1.0 w := math.Abs(x2 - x1) lastBisectWidth := w var nFalseP int var x3, f3, bestX, bestF float64 for i := 0; i < maxIterations; i++ { switch state { case bisect: x3 = 0.5 * (x1 + x2) if x3 == x1 || x3 == x2 { // i.e., x1 and x2 are successive floating-point numbers bestX = x3 if x3 == x1 { bestF = f1 } else { bestF = f2 } return fSolveConverged, bestX, bestF, w } f3 = f(x3, fExtra) if f3 == 0 { return fSolveExact, x3, f3, w } if f3*f2 < 0 { x1 = x2 f1 = f2 } x2 = x3 f2 = f3 w = math.Abs(x2 - x1) lastBisectWidth = w if bisectTil > 0 { if w < bisectTil { bisectTil = -1.0 gamma = 1.0 nFalseP = 0 state = falseP } } else { gamma = 1.0 nFalseP = 0 state = falseP } case falseP: s12 := (f2 - gamma*f1) / (x2 - x1) x3 = x2 - f2/s12 f3 = f(x3, fExtra) if f3 == 0 { return fSolveExact, x3, f3, w } nFalseP++ if f3*f2 < 0 { gamma = 1.0 x1 = x2 f1 = f2 } else { // Anderson-Bjorck method g := 1.0 - f3/f2 if g <= 0 { g = 0.5 } gamma *= g } x2 = x3 f2 = f3 w = math.Abs(x2 - x1) // Sanity check. For every 4 false position checks, see if we really are // decreasing the interval by comparing to what bisection would have // achieved (or, rather, a bit more lenient than that -- interval // decreased by 4 instead of by 16, as the fp could be decreasing gamma // for a bit). Note that this should guarantee convergence, as it makes // sure that we always end up decreasing the interval width with a // bisection. if nFalseP > bisectIter { if w*bisectWidth > lastBisectWidth { state = bisect } nFalseP = 0 lastBisectWidth = w } } tol := absErr + relErr*math.Max(math.Max(math.Abs(x1), math.Abs(x2)), 1.0) if w <= tol { if math.Abs(f1) < math.Abs(f2) { bestX = x1 bestF = f1 } else { bestX = x2 bestF = f2 } return fSolveConverged, bestX, bestF, w } } return fSolveMaxIterations, x3, f3, w }