Various numerical methods.
For demos of some of these algorithms refer to the following notebook.
function bisect (f:string, a:number, b:number, y:number=0) Uses the bisection method to solve f(x)=y in [a,b] f is either a function of one real variable, or an expression in x. Addtional parameters to "bisect" after a semicolon will be passed to the function f. The bisection routine assumes a zero in the interval, i.e., it assumes different signs of f(a) and f(b). The routine bisects the interval until the desired accuracy is reached. It is a very stable function. It is possible to specify the accuracy epsilon with eps=... as last parameter. f : function or expression in x a,b : interval bounds y : target value >bisect("x^2",0,2,y=2) 1.41421356237 See:
newton (Numerical Algorithms),
newton (Maxima Documentation),
secant (Numerical Algorithms),
ibisect (Interval Solvers and Guaranteed Solutions)
function secant (f:string, a:complex scalar, b=none, y:number=0) Uses the secant method to solve f(x)=y in [a,b] f is either a function of one real variable, or an expression in x. Additional parameters to "secant" after a semicolon will be passed to the function f. The secant routine uses an approximation of the derivative to find a better approximation to the solution. This function is almost as good as the Newton method, but does not require the derivative of f. Note that the result is not guaranteed to be within the interval [a,b]. The function may fail, if the approximated derivative is too close to 0. The parameter b is optional. In this case, the function starts with an interval around a. It is possible to specify the accuracy epsilon with eps=... as last parameter. f : function or expression in x a,b : interval bounds, b is optional y : target value >bisect("x^2",0,2,y=2) 1.41421356237 > >secant("x^2",1,y=2) 1.41421356237 >function f(x,a) := x^2-a >secant("f",5;3), %^2 1.73205080757 3 See:
bisect (Numerical Algorithms),
newton (Numerical Algorithms),
newton (Maxima Documentation)
function solve (f:string, a:complex scalar, b=none, y:number=0) Solve f(x)=y near a. Alias to "secant". See:
secant (Numerical Algorithms)
function root (f:string, %x) Find the root of an expression f by changing the variable x. This is a function for interactive use. f should be an expression in several variables. All variables but x are constant, and only x is changed to solve f(x,...)=0. Other variables must be set globally. Note that the actual name of the variable x may be different from "x", so you can solve for any variable in the expression f. All variables in the expression must be global! f : expression in several variables %x : variable to bo solved for >a=2; x=2; root("x^x-a",x); "a = "+a, "x = "+x, a = 2 x = 1.55961046946 >a=2; x=2; root("x^x-a",a); "a = "+a, "x = "+x, a = 4 x = 2 See:
secant (Numerical Algorithms)
function simpson (f:string, a:number, b:number, n:natural=50, maps:integer=0) Integrates f in [a,b] using the Simpson method The Simpson method is stable, but not fast. The results are not accurate. The method is exact for polynomials of degree 3. Use "gauss" for better results. f is either a function of one real variable, or an expression in x. Addtional parameters to "bisect" after a semicolon will be passed to the function f. f : function or expression in x a,b : interval bounds n : number of points to use maps : flag, if the expression needs to be mapped (0 or 1). Functions are always mapped. See:
romberg (Numerical Algorithms),
romberg (Maxima Documentation),
gauss (Gauss Qudrature)
function simpsonfactors (n) Returns the vector 1,4,2,...,2,4,1.
function romberg (f:string, a:number, b:number, m:natural=10, maps:integer=0) Romberg integral of f in [a,b] This method gets very accurate results for many functions. The Romberg algorithm uses an extrapolation of the trapezoidal rule for decreasing step sizes. However, the Gauss method is to be preferred for functions, which are close to polynomials. The function may fail, if f is not smooth enough. f is either a function of one real variable, or an expression in x. Addtional parameters to "bisect" after a semicolon will be passed to the function f. f : function or expression in x a,b : interval bounds m : number of points to use for the start of the algorithm maps : flag, if the experssion needs to be mapped (0 or 1). Functions are always mapped. See:
simpson (Numerical Algorithms),
gauss (Gauss Qudrature),
integrate (Differential Equations),
integrate (Maxima Documentation)
function iterate (f:string, x:numerical, n:natural=0) Iterates f starting from x to convergence. This function iterates xnew=f(x), until xnew~=x, or the maximal number of iterations n is exceeded. It returns the limit for n=0 (the default), or n column vectors. f is either a function of a scalar returning a scalar, or a function of a vector returning a vector. The function will work with column or row vectors, and even for matrices, if n==0. f can also be an expression of a variable x. f : function or expression x : start point for the algorithm n : optional number of steps
function niterate (f:string, x:numerical, n:natural) Iterates f n times starting from x. Works like iterate with n>0. f : function or expression x : start point for the algorithm n : number of steps
function sequence (f:string, x:numerical, n:natural) Computes sequences recursively. This functions is used to compute x[n]=f(x[1],...,x[n-1]) recursively. It will work for scalar x[i], or for column vectors x[i]. Depending on the iteration function, more than one start value x[1] to x[k] is needed. These start values are stored in a row vector x0, or for vector iteration in a matrix x0 as columns. The iteration count n determines the highest x[n], which is computed. f is a function f(x,n), or an expression in x an n. The parameter x contains the values x[1],...,x[n-1] computed so far as elements of a row vector, or as columns in a matrix x. In case of a vector sequence, the function can also work with row vectors, depending on the result of the function f. Note, that the start value must fit to the result of the function. f : function f(x,n) (x : mx(n-1) or (n-1)xm, result: mx1 or 1xm) x : start values (mxk or kxm) n : desired number of values See:
iterate (Numerical Algorithms)
function newton (f:string, df:string, x:complex scalar, y:number=0) Solves f(x)=y with the Newton method. The Newton method needs the derivative of the function f. This derivative must be computed by the function df. The Newton method will start in x, and stop if the desired accuracy is reached. It is possible to set this accuracy with eps=... The function will only work for real or complex variables. For systems of equations, see "newton2". For intervals, see "inewton". f and df must be a function of a real or complex scalar, returning a scalar, or an expression of x. f : function or expression in x (real or complex, scalar) df : function or expression in x (real or complex, scalar) x : start value y : target value See:
secant (Numerical Algorithms),
bisect (Numerical Algorithms),
newton2 (Numerical Algorithms),
inewton (Interval Solvers and Guaranteed Solutions),
mxminewton (Maxima Functions for Euler)
function newton2 (f:string, df:string, x:numerical) Newton methods to solve a system of equations. The mulitdimensional Newton method can solve the equation f(v)=0 for row vectors v. Here, f must map row vectors to row vectors. This function needs the Jacobian matrix of f, which is provided throught the function df. The Newton method will fail, if the Jacobian gets singular during the computation. f must be a function of a row vector, returning a row vector, or an expression of the row vector x. df must be a function of a row vector, returning the Jacobian matrix of f. Additional parameters after a semicolon are passed to the functions f and df. The function f can also use column vectors. In this case, the start value must be a column vector. f : function f(v) (v: 1xn or nx1, result 1xn or nx1) df : function df(v) (v: 1xn or nx1, result: nxn) x : start point (1xn or nx1) See:
broyden (Numerical Algorithms),
neldermin (Numerical Algorithms),
inewton2 (Interval Solvers and Guaranteed Solutions)
function steffenson (f:string, x:scalar, n:natural=0) Does n steps of the Steffenson operator, starting from x. This is a similar function as iterate. However, the function assumes an asymptotic expansion of the convergence, and uses the Steffenson operator to speed up the convergence. f is either a function of a scalar returning a scalar, or a function of a column vector returning a column vector. f can also be an expression of a variable x. f : function or expression of a scalar x x : start point n : optional number of iterations See:
iterate (Numerical Algorithms)
function nsteffenson (ggg:string, x0:scalar, n:natural) Does n steps of the Steffenson operator, starting from x0. Works like "steffenson", but returns the history of the iteration. See:
niterate (Numerical Algorithms)
function aitkens (x: vector) Improves the convergence of the sequence x. The Aitkens operator assumes that a-x[n] has a geometric convergence 1/c^n to 0. With this information, it is easy to compute the limit a. x : row vector (1xn, n>2) Return a row vector (1x(n-2)) with faster convergence. See:
steffenson (Numerical Algorithms)
DifMatrix:=%setupdif(5);
function diff (f:string, x:numerical, n:natural=1, h:number=epsilon()^(1/3)) Compute the n-th (n<6) derivative of f in the points x. Numerical differentiation is somewhat inaccurate. To get a good approximation, the first derivative uses 4 function evaluations. f is either a function in a scalar, or an expression in x. Additional parameters after a semicolon are passed to f. f : function or expression a scalar x x : the point, where the derivative should be computed n : degree of derivative h : optional delta
function dif (f, x, n=1, h=epsilon()^(1/3)) Compute the n-th (n<6) derivative of f in the points x. Alias to diff.
function fftfold (v:vector, w:vector) fold v with w using fast Fourier transformation. The length of the vector v should have small prime factors. For large v, this is much faster than fold(v,w). However, the vector v is folded as a periodic vector. To get the same result as fold(v,w), delete the first cols(w)-1 coefficients. v : row vector (1xn) w : row vector (1xm, m<=n) See:
fold (Euler Core)
The Broyden algorithm is a Quasi-Newton method for solving F(x)=0 with vectors x. It does not need a derivative. Instead, it uses an approximation of the Jacobian Df to start with, and updates it at each step.
function broyden (f:string, xstart:real, A:real=none) Finds a zero of f with the Broyden algorithm. The Broyden algorithm is an iterative algorithm just like the Newton method, which solves the equation f(v)=0 for vectors v. It tries to approximate the Jacobian of the function f using the previous iteration steps. The algorithm will fail, if the Jacobian of f in the zero is singular. The function f must take a vector v, and return a vector. Additional parameters to "broyden" after the semicolon are passed to f. The function can work with column or row vectors. The start vector must be of the same form. x is the start value, and the optional matrix A is a approximation of the Jacobian matrix of f. To change the accuracy, you can specify an optional eps=... as last parameter. f : function of one vector xstart : start point, a row vector A : optional start for the Jacobian Returns the solution as a vector. See:
nbroyden (Numerical Algorithms)
function nbroyden (f:string, xstart:real vector, nmax:natural, A:real=none) Do nmax steps of the Broyden algorithm. This function works like "broyden", but returns the steps of the algorithm. Moreover, it does at most nmax steps. f : function of one row vector xstart : start point, a row vector nmax : maximal number of steps A : optional start for the Jacobian Returns a matrix, with one step of the algorithm in each row, and the current approximation of the Jacobian.
function hermiteinterp (x:vector, y) Divided differences for Hermite interpolation. x : row vector of interpolation points. x can contain double points in immediate succession. These points use derivatives for the interpolation. y : Either a function df(x,n), which computes the n-th derivatives of f (including the 0-th derivative), or a vector of values and derivatives in the form f(x),f'(x),f''(x),... for a multiple interpolation point x. Returns the divided differences, which can be used in interpval (divdifeval) as usual. The multiplicities in x can be generated with multdup. >xp=multdup(1:3,2) [ 1 1 2 2 3 3 ] >yp=[1,0,0,0,-1,0] [ 1 0 0 0 -1 0 ] >d=hermiteinterp(xp,yp) [ 1 0 -1 2 -1.5 1.5 ] >plot2d("interpval(xp,d,x)",0.9,3.1); >function fh (x,n,f,fd) ... $ if n==0 then return f(x); $ elseif n==1 then return fd(x); $ else error("n=0 or n=1!"); $ endfunction >xp=multdup(chebzeros(0,1,5),2); >expr &= sqrt(x); >d=hermiteinterp(xp,"fh";expr,&diff(expr,x)); >plot2d("sqrt(x)-interpval(xp,d,x)",0,1); See:
divdif (Euler Core),
divdifeval (Euler Core),
multdup (Numerical Algorithms)
function hermitedivdif (x:vector, y) Divided differences for Hermite interpolation. Alias to hermiteinterp See:
hermiteinterp (Numerical Algorithms)
function multdup (x:numerical, n:nonnegative integer) Repeat the rows or columns of x with multiplicities in n E.g., if n[1]=2, and n is a row vector, the first column of x is duplicated. If n is a column vector, the rows of x are duplicated, if it is a row vector, the columns are duplicated. n can be shorter than the number of rows resp. columns of x. If n is scalar, it acts on the rows, duplicating the first row n times. See:
dup (Euler Core),
hermite (Maxima Documentation)
function remez (x:vector, y:vector, deg:index, tracing:integer=0, remezeps=0.00001) Best polynomial approximation. The Remez algorithm computes the best polynomial approximation to the data (x[i],y[1]) of degree deg. This algorithm uses a simultanous exchange, which requires the points x[i] to be sorted. x : vector of points on the real line, sorted. y : vector of values at these points. deg : degree of the polynomial. tracing : if non-zero, the steps will be plotted. Returns {t,d,h,r} t : the alternation points d : the divided difference form of the best approximant h : the discrete error (with sign) r : the rightmost alternation point, which is missing in t. To evaluate the polynomial in v use interpval(t,d,v).
function fmin (f:string, a:real scalar, b:real scalar) Compute the minimum of the convex function f in [a,b]. Uses the golden cut method, starting with the interval [a,b]. The method takes about 60*(b-a) function evalations for full accuracy. f is either an expression in x, or a function of x. Additional parameters are passed to a function f. You can specify an epsilon eps with eps=... as last parameter. See:
fmax (Numerical Algorithms)
function fmax (f:string, a:real scalar, b:real scalar) Compute the maximum of the concave function f in [a,b]. Uses the golden cut method, starting with the interval [a,b]. The method takes about 60*(b-a) function evalations for full accuracy. f is either an expression in x, or a function of x. Additional parameters after a semicolon are passed to a function f. You can specify an epsilon eps with eps=... as last parameter. See:
fmin (Numerical Algorithms)
function fextrema (f:string, a:real scalar, b:real scalar, n:integer=100) Find all internal extrema of f in [a,b]. f may be an expression in x or a function. Additional parameters after a semicolon are passed to a function f. You can specify an epsilon eps with eps=... as last parameter. Returns {minima,maxima} (vectors, possibly of length 0) See:
fmax (Numerical Algorithms),
fmin (Numerical Algorithms)
function brentmin (f:string, a:real scalar, d=0.1, eps=epsilon()) Compute the minimum of f around a with Brent's method. d is an initial step size to look for a starting interval. eps is the final accuracy. Returns the point of minimal value. f is an expression in x, or a function in f. Additional parameters after a semicolon are passed to a function f. Return the point of minimum. See:
fmin (Numerical Algorithms),
fmax (Numerical Algorithms)
function neldermin (f:string, v:real, d=0.1, eps=epsilon(), tries=2) Compute the minimum of f around v with the Nelder-Mead method. The Nelder-Mead method is a stable, but slow method to search for a local minimum of a function without using any derivative. It should not be used for high dimensions. Of course, it can also be applied to solve a system equations by minimizing the norm of the errors. f must be function of a row or column vector x, returning a real value. Addditional parameters after the semicolon are passed to f. f can also be an expression depending on a vector x. d is an optional initial step size for a starting simplex. eps is the final accuracy. f : function f(v) (v : 1xn or nx1, result: scalar) v : start point for the search (1xn or nx1) d : optional size of start simplex eps : optional accuracy tries : number of restarts of the algorithm Return the point of minimum (1xn vector). See:
nelder (Euler Core)
function nelder (f:string, v:real, d=0.1, eps=epsilon()) Compute the minimum of f around v with the Nelder-Mead method. Alias to neldermin. See:
neldermin (Numerical Algorithms)
function nlfit (f,Df,v) Minimizes f(x) from start vector v. This method is named after Levenberg-Marquardt. f(x) maps 1xn vectors to 1xm vectors (m>n) Df(x) is the Jacobian of f.
function newtonbarrier (f:string, df:string, Hf:string, A:real, b:real column, x:real vector, lambda:positive real=1, c:positive real=0.1, cfactor:positive real=0.1, history=false, eps=epsilon()) Newton-Barrier method starting from interior x. This function can minimize a convex function function f(x) subject to conditions Ax <= b. It needs the gradient df and the Hessian Hf of f as functions. The start point x must be an interior point with Ax < b. Note: f must be convex. f,df,Hf : Functions taking a vector as a parameter. A,b : Conditions Ax<=b x : Start point with Ax<b lambda,c,cfactor : Variables for the algorithms. You can take the default. history : If true, the function returns a matrix with one step x in each row. Else it returns the last x only. eps : The absolute accuracy of the solution. The algorithm returns, when |x-xnew|<epsilon.
For an example of the Newton Barrier Method see the following introduction notebook.
function spline (x,y) Defines the natural Spline at points x(i) with values y(i). The natural spline is the spline using cubic polynomials between the points x[i], smooth second derivative, and linear outside the interval x[1] to x[n]. The points x[i] must be sorted. The function returns the second derivatives at these points. With this information, the spline can be evaluated using "splineval". >xp=1:10; yp=intrandom(1,10,3); s=spline(xp,yp); >plot2d("splineval(x,xp,yp,s)",0,11); >plot2d(xp,yp,>points,>add); See:
splineval (Numerical Algorithms)
function splineval (t:number; x:vector, y:vector, h:vector) Evaluates the cubic interpolation spline for (x(i),y(i)) with second derivatives h(i) at t
function evalspline (t:number, x:vector, y:vector, s:vector) Evaluates the cubic interpolation spline for (x(i),y(i)) with second derivatives s(i) at t See:
splineval (Numerical Algorithms)
function linsplineval (x:number; xp:vector, yp:vector, constant:integer=1) Evaluates the linear interpolating spline for (xp[i],yp[i]) at x xp must be sorted. Outside of the points of xp, the spline is continued as in the closest interval to x, or as a constant function. >xp=1:10; yp=intrandom(1,10,3); >plot2d("linsplineval(x,xp,yp,<constant)",0,11); >plot2d(xp,yp,>points,>add);
function evallinspline (x:number, xp:vector, yp:vector, constant:integer=1) Evaluates the linear interpolating spline for (xp[i],yp[i]) at x See:
linsplineval (Numerical Algorithms)