blob: 490ca8ac7e0c0847f573769e34d12033b9ff4d95 [file] [log] [blame] [edit]
// 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
}