iconEuler Reference

Numerical Algorithms

Various numerical methods.

For demos of some of these algorithms refer to the following notebook.

Solve Equations

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)

Integrals

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)

Iterations

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)

Newton

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)

Differentiation

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.

Folding

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)

Broyden Algorithm

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.
  

Hermite Interpolation

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)

Remez Algorithm

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).

Nonlinear Optimization

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.

Splines

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)

Documentation Homepage