# README
Gosl. num. Fundamental numerical methods
This package implements essential numerical methods such as for root finding, numerical quadrature, numerical differentiation, and solution of simple nonlinear problems.
While the sub-package num/qpck provides advanced quadrature schemes (by wrapping Quadpack), this package implements few (simpler) methods to compute numerical integrals. Here, there are two kinds of algorithms: (1) basic methods for discrete data; and (2) using refinement for integrating general functions.
Example: Using Brent's method:
Find the root of
y(x) = x³ - 0.165 x² + 3.993e-4
within [0, 0.11]. See figure below. Note: we have to make sure that the root is bounded otherwise Brent's method doesn't work.
// y(x) function
yx := func(x float64) float64 {
return math.Pow(x, 3.0) - 0.165*math.Pow(x, 2.0) + 3.993e-4
}
// range: be sure to enclose root
xa, xb := 0.0, 0.11
// initialize solver
solver := num.NewBrent(yx, nil)
// solve
xo := solver.Root(xa, xb)
// output
yo := yx(xo)
io.Pf("x = %v\n", xo)
io.Pf("f(x) = %v\n", yo)
io.Pf("nfeval = %v\n", solver.NumFeval)
io.Pf("niter. = %v\n", solver.NumIter)
Output of Brent's solution:
it x f(x) err
1.0e-14
0 1.100000000000000e-01 -2.662000000000001e-04 5.500000000000000e-02
1 6.600000000000000e-02 -3.194400000000011e-05 3.300000000000000e-02
2 6.044444444444443e-02 1.730305075445823e-05 2.777777777777785e-03
3 6.239640011030302e-02 -1.676981032316081e-07 9.759778329292944e-04
4 6.237766369176578e-02 -7.323468182796403e-10 9.666096236606754e-04
5 6.237758151338346e-02 3.262039076357137e-15 4.108919116063703e-08
6 6.237758151374950e-02 0.000000000000000e+00 4.108900814037142e-08
x = 0.0623775815137495
f(x) = 0
nfeval = 8
niter. = 6
Example: Using Newton's method:
Same problem as before.
// Function: y(x) = fx[0] with x = xvec[0]
fcn := func(fx, xvec la.Vector) {
x := xvec[0]
fx[0] = math.Pow(x, 3.0) - 0.165*math.Pow(x, 2.0) + 3.993e-4
}
// Jacobian: dfdx(x) function
Jfcn := func(dfdx *la.Matrix, x la.Vector) {
dfdx.Set(0, 0, 3.0*x[0]*x[0]-2.0*0.165*x[0])
return
}
// trial solution
xguess := 0.03
// initialize solver
neq := 1 // number of equations
useDn := true // use dense Jacobian
numJ := false // numerical Jacobian
var o num.NlSolver
o.Init(neq, fcn, nil, Jfcn, useDn, numJ, nil)
// solve
xvec := []float64{xguess}
o.Solve(xvec, false)
// output
fx := []float64{123}
fcn(fx, xvec)
xo, yo := xvec[0], fx[0]
io.Pf("x = %v\n", xo)
io.Pf("f(x) = %v\n", yo)
io.Pf("nfeval = %v\n", o.NFeval)
io.Pf("niter. = %v\n", o.It)
Output of NlSolver:
it Ldx fx_max
(1.0e-04) (1.0e-09)
0 0.000000000000000e+00 2.778000000000000e-04
1 3.745954692556634e+06 5.421253067129628e-05
2 6.176571448942142e+05 1.391803634400563e-06
2 1.515117884960284e+04 5.314115983194589e-10
. . . converged with fx_max. nit=2, nFeval=4, nJeval=3
x = 0.062377521883073835
f(x) = 5.314115983194589e-10
nfeval = 4
niter. = 2
Example: Quadrature with discrete data
Example: Quadrature with methods using refinement
Example: numerical differentiation
Check first and second derivative of y(x) = sin(x)
// define function and derivative function
yFcn := func(x float64) float64 { return math.Sin(x) }
dydxFcn := func(x float64) float64 { return math.Cos(x) }
d2ydx2Fcn := func(x float64) float64 { return -math.Sin(x) }
// run test for 11 points
X := utl.LinSpace(0, 2*math.Pi, 11)
io.Pf(" %8s %23s %23s %23s\n", "x", "analytical", "numerical", "error")
for _, x := range X {
// analytical derivatives
dydxAna := dydxFcn(x)
d2ydx2Ana := d2ydx2Fcn(x)
// numerical derivative: dydx
dydxNum := num.DerivCen5(x, 1e-3, func(t float64) float64 {
return yFcn(t)
})
// numerical derivative d2ydx2
d2ydx2Num := num.DerivCen5(x, 1e-3, func(t float64) float64 {
return dydxFcn(t)
})
// check
chk.PrintAnaNum(io.Sf("dy/dx @ %.6f", x), 1e-10, dydxAna, dydxNum, true)
chk.PrintAnaNum(io.Sf("d²y/dx² @ %.6f", x), 1e-10, d2ydx2Ana, d2ydx2Num, true)
}
Examples: nonlinear problems
Find x0 and x1 such that f0 and f1 are zero, with:
f0(x0,x1) = 2.0*x0 - x1 - exp(-x0)
f1(x0,x1) = -x0 + 2.0*x1 - exp(-x1)
Using analytical (sparse) Jacobian matrix
Output:
it Ldx fx_max
(1.0e-05) (1.0e-15)
0 0.000000000000000e+00 4.993262053000914e+00
1 8.266404824090484e+09 9.204814001140181e-01
2 7.824673760247719e+08 9.107574803964624e-02
3 9.482829646746747e+07 9.541527544986161e-04
4 1.014523823737919e+06 1.051153347697564e-07
5 1.117908116077260e+02 1.221245327087672e-15
5 1.298802024636321e-06 1.110223024625157e-16
. . . converged with fx_max. nit=5, nFeval=12, nJeval=6
x = [0.5671432904097838 0.5671432904097838] expected = [0.5671 0.5671]
f(x) = [-1.1102230246251565e-16 -1.1102230246251565e-16]
Using numerical Jacobian matrix
Output:
it Ldx fx_max
(1.0e-05) (1.0e-15)
0 0.000000000000000e+00 4.993262053000914e+00
1 8.266404846831476e+09 9.204814268661379e-01
2 7.824673975180904e+08 9.107574923794759e-02
3 9.482829862677714e+07 9.541518971115659e-04
4 1.014522914054942e+06 1.051134990159852e-07
5 1.117888589104628e+02 1.554312234475219e-15
5 1.653020764599351e-06 1.110223024625157e-16
. . . converged with fx_max. nit=5, nFeval=24, nJeval=6
xx = [0.5671432904097838 0.5671432904097838] expected = [0.5671 0.5671]
f(xx) = [-1.1102230246251565e-16 -1.1102230246251565e-16]
Using analytical dense Jacobian matrix
Just replace Jfcn with
JfcnD := func(dfdx [][]float64, x []float64) error {
dfdx[0][0] = 2.0+math.Exp(-x[0])
dfdx[0][1] = -1.0
dfdx[1][0] = -1.0
dfdx[1][1] = 2.0+math.Exp(-x[1])
return nil
}
API
# Functions
CompareJac compares Jacobian matrix (e.g.
CompareJacDense compares Jacobian matrix (e.g.
DerivBwd4 approximates the derivative df/dx using backward differences with 4 points.
DerivCen5 approximates the derivative df/dx using central differences with 5 points.
DerivFwd4 approximates the derivative df/dx using forward differences with 4 points.
EqCubicSolveReal solves a cubic equation, ignoring the complex answers.
GaussJacobiXW computes positions (xi) and weights (wi) to perform Gauss-Jacobi integrations.
GaussLegendreXW computes positions (xi) and weights (wi) to perform Gauss-Legendre integrations Input: x1 -- lower limit of integration x2 -- upper limit of integration n -- number of points for quadrature formula Reference: [1] Press WH, Teukolsky SA, Vetterling WT, Fnannery BP (2007) Numerical Recipes: The Art of Scientific Computing.
Jacobian computes Jacobian (sparse) matrix Calculates (with N=n-1): df0dx0, df0dx1, df0dx2, ..
LineSearch finds a new point x along the direction dx, from x0, where the function has decreased sufficiently.
LinFit computes linear fitting parameters.
LinFitSigma computes linear fitting parameters and variances (σa,σb) in the estimates of a and b Errors on y-direction only
y(x) = a + b⋅x
See page 780 of [1] Reference: [1] Press WH, Teukolsky SA, Vetterling WT, Fnannery BP (2007) Numerical Recipes: The Art of Scientific Computing.
NewBracket returns a new bracket-er object.
NewBrent returns a new Brent structure ffcn -- function f(x) Jfcn -- derivative df(x)/dx [optional / may be nil].
NewLineSolver returns a new LineSolver object size -- length(x) ffcn -- scalar function of vector: y = f({x}) Jfcn -- vector function of vector: {J} = df/d{x} @ {x} [optional / may be nil].
NewNlSolver creates a new NlSolver F is the f(x) function f:vector, x:vector Will use numerical Jacobian (with sparse solver) by default.
NewNlSolverConfig creates a new NlSolverConfig Default values: CteJac = false LinSearch = false LinSchMaxIt = 20 MaxIt = 20 ChkConv = false Atol = 1e-8 Rtol = 1e-8 Ftol = 1e-9.
QuadCs performs automatic integration (quadrature) using the cosine or sine weights QUADPACK routine AWOE (Automatic with weight, Oscillatory)
INPUT: a -- lower limit of integration b -- upper limit of integration ω -- omega useSin -- use sin(ω⋅x) instead of cos(ω⋅x) fid -- index of goroutine (to avoid race problems) f -- function defining the integrand
OUTPUT: b b res = ∫ f(x) ⋅ cos(ω⋅x) dx or res = ∫ f(x) ⋅ sin(ω⋅x) dx a a
.
QuadDiscreteSimps2d approximates a double integral over the x-y plane with the elevation given by data points f[npts][npts].
QuadDiscreteSimpsonRF approximates the area below the discrete curve defined by [xa,xy] range and y function.
QuadDiscreteTrapz2d approximates a double integral over the x-y plane with the elevation given by data points f[npts][npts].
QuadDiscreteTrapzRF approximates the area below the discrete curve defined by [xa,xy] range and y function.
QuadDiscreteTrapzXF approximates the area below the discrete curve defined by x points and y function.
QuadDiscreteTrapzXY approximates the area below the discrete curve defined by x and y points.
QuadExpIx approximates the integral of f(x) ⋅ exp(i⋅m⋅x) with i = √-1
INPUT: a -- lower limit of integration b -- upper limit of integration m -- coefficient of x fid -- index of goroutine (to avoid race problems) f -- function defining the integrand
OUTPUT: b b b res = ∫ f(x) ⋅ exp(i⋅m⋅x) dx = ∫ f(x) ⋅ cos(m⋅x) dx + i ⋅ ∫ f(x) ⋅ sin(m⋅x) dx a a a
.
QuadGaussL10 approximates the integral of the function f(x) between a and b, by ten-point Gauss-Legendre integration.
QuadGen performs automatic integration (quadrature) using the general-purpose QUADPACK routine AGSE (Automatic, general-purpose, end-points singularities).
SecondDerivCen3 approximates the second derivative d²f/dx² using central differences with 3 points.
SecondDerivCen5 approximates the second derivative d²f/dx² using central differences with 5 points.
SecondDerivMixedO2 approximates ∂²f/∂x∂y @ x={x,y} using O(h²) formula.
SecondDerivMixedO4v1 approximates ∂²f/∂x∂y @ x={x,y} using O(h⁴) formula from http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/central-differences.
SecondDerivMixedO4v2 approximates ∂²f/∂x∂y @ x={x,y} using O(h⁴) formula from http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/central-differences.
# Variables
smallest number satisfying 1 + EPS > 1.
# Structs
Bracket implements routines to bracket roots or optima
A root of a function is known to be bracketed by a pair of points, a and b, when the function has opposite sign at those two points.
Brent implements Brent's method for finding the roots of an equation.
ElementarySimpson structure implements the Simpson's method for quadrature with refinement.
ElementaryTrapz structure is used for the trapezoidal integration rule with refinement.
LineSolver finds the scalar λ that zeroes or minimizes f(x+λ⋅n).
NlSolver implements a solver to nonlinear systems of equations References: [1] G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical computations.
NlSolverConfig holds the configuration input for NlSolver.
# Interfaces
QuadElementary defines the interface for elementary quadrature algorithms with refinement.