fuchsia / third_party / github.com / gonum / gonum / refs/heads/upstream/asm/copy / . / mathext / roots.go

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

} |