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