iconEuler Reference

Numerical Algorithms

Content

Various numerical methods.

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

Solving Numerical Equations

function bisect (f$:call, a:number, b:number, y:number=0,
    eps=epsilon())
  Bisection 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 "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 an
  additional parameter. The default is the internal epsilon. The
  absolute error is used.
  
  There is the interval alternative ibisect().
  
  f : function or expression in x
  a,b : interval bounds
  y : target value
  eps : epsilon for the absolute error.
  
  >bisect("x^2",0,2,y=2) // solve x^2=2
  1.41421356237
  >bisect("x^2",0,2,y=2,eps=0.001)
  1.41455078125
  
  >function f(x,a) := exp(-a*x)*cos(x)
  >x=bisect("f",0,pi/2;2,y=0.5) // solve f(x,2)=0.5
  0.320449786201
  >f(x,2) // test
  0.5
  
  See: 
newton (Numerical Algorithms),
newton (Maxima Documentation),
secant (Numerical Algorithms),
ibisect (Numerical Algorithms),
call (Euler Core)
function intbisect (f$:call, a:number, b:number, y:number=0)
  Bisection method to find an integer solution for f(x)=y.
  
  This function should only be used for monotone functions. It makes
  sense for functions, which are defined only for integer values.
  
  The solution satisfies f(floor(x))>=y, f(ceil(x))<=y (reversed if
  f is increasing). In the last step, x interpolates f in a linear
  way between floor(x) and ceil(x).
  
  y : target value
  
  >intbisect("chidis(6.5,x)",1,10,y=0.5)
  7.15896300259
  >chidis(6.5,[7,8])
  [0.517277,  0.408592]
  
  See: 
call (Euler Core)
function secant (f$:call, a:complex scalar, b=none,
    y:number=0, eps=none)
  Secant method to solve f(x)=y in [a,b]
  
  f is either a function of one real variable, or an expression or a
  call collection 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]. For an alternative, see secantin(). Moreover, 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.
  
  The method can also be used for complex functions.
  
  It is possible to specify the accuracy 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) // expression
  1.41421356237
  >secant({{"x^2-a",a=2}},1) // call collection
  1.41421356237
  >function f(x,a) := x^2-a
  >secant("f",5;3), %^2
  1.73205080757
  3
  >secant("x^2",1,y=I) // complex equation x^2=I
  0.707107+0.707107i
  
  See: 
bisect (Numerical Algorithms),
newton (Numerical Algorithms),
newton (Maxima Documentation),
secantin (Numerical Algorithms),
call (Euler Core)
function secantin (f$:call, a:real scalar, b:real scalar,
    y:number=0,eps=epsilon())
  Secant method to solve f(x)=y in [a,b]
  
  Works like the secant method. But the result is guaranteed to be
  in [a,b]. Use this for functions, which are not defined outside of
  [a.b]. The parameter "eps" is used for the absolute error.
  
  >secant("sqrt(x)*exp(-x)",0.4,y=0.3) // first leads to x<0
  Floating point error!
  Error in sqrt
  Error in expression: sqrt(x)*exp(-x)
  secant:         x0=x1; y0=y1; x1=x2; y1=f$(x2,args())-y;
  >secantin("sqrt(x)*exp(-x)",0,0.4,y=0.3) // solve in [0,0.4]
  0.112769901579
  
  See: 
bisect (Numerical Algorithms),
newton (Numerical Algorithms),
newton (Maxima Documentation),
call (Euler Core)
function solve (f$:call, a:complex scalar, b=none,
    y:number=0, eps=none)
  Solve f(x)=y near a.
  
  Alias to "secant".
  
  See: 
secant (Numerical Algorithms),
call (Euler Core)
function root (f$:string, %x, eps=none)
  Find the root of an expression f by changing the variable x.
  
  This is a function for interactive use. It will not work properly
  inside of functions. Use solve() instead.
  
  f must be an expression in several variables. All variables but x
  are constant, and only x is changed to solve f(x,...)=0. All
  variables in the expression must be global variables. 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. But the parameter must
  be a variable.
  
  f : expression in several variables
  %x : variable to be 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),
solve (Numerical Algorithms),
solve (Maxima Documentation),
newton (Numerical Algorithms),
newton (Maxima Documentation)

Numerical Iteration

function iterate (f$:call, x:numerical, n:natural=0,
    till=none, eps=none, best=false)
  Iterates f starting from x to convergence.
  
  This function iterates xnew=f(x), until xnew~=x (using the
  internal epsilon or the optional eps), or the maximal number of
  iterations n is exceeded. If >best is set the iteration will
  continue until abs(x-f(x)) does no longer get smaller. It returns
  the limit for n=0 (the default), or n iterated values. There is a
  "till" expression, which can also stop the iteration.
  
  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 or a call collection of a variable x.
  
  The iteration can also be used for intervals or interval vectors.
  In this case, The iteration will stop until the left and right
  interval bounds are closer than epsilon. If >best is set the
  iteration continues until the interval length does not get any
  smaller.
  
  f : function or expression
  x : start point for the algorithm
  n : optional number of steps (0 = till convergence)
  till : optional end condition
  best : Iterate until best result is achieved (takes long if the
  limit is 0).
  
  >iterate("cos",1) // solves cos(x)=x, iterating from 1
  0.739085133216
  >longest iterate("(x+2/x)/2",1,n=7)'
  1
  1.5
  1.416666666666667
  1.41421568627451
  1.41421356237469
  1.414213562373095
  1.414213562373095
  1.414213562373095
  >iterate("x+1",1000,till="isprime(x)")
  1009
  >function f(x,A,b) := A.x-b
  >A=random(3,3)/2; b=sum(A);
  >iterate("f",zeros(3,1);A,b)
  -11.7766
  -12.4523
  -14.938
  >iterate("(middle(x)+2/x)/2",~1,2~) // interval iteration
  ~1.414213562372,1.414213562374~
  >iterate("(x+I/x)/2",1) // complex iteration
  0.707107+0.707107i
  
  Returns {x,n}, where n is the number of iterations needed
  
  See: 
niterate (Numerical Algorithms),
call (Euler Core)
function niterate (f$:call, x:numerical, n:natural, till=none)
  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 iterations
  till : optional end condition
  
  Returns v, where v is the vector of iterations and n the
  number of iterations.
  
  See: 
call (Euler Core)
function sequence (f$:call, x:numerical, n:natural)
  Computes sequences recursively.
  
  This function can be 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 in a matrix for vector iteration). The iteration
  counter n can be used in the expression.
  
  f$ must be a function f$(x,n), or an expression or a call
  collection 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) where x is m x (n-1) or (n-1) x m,
  and the result m x 1 or 1 x m
  x : start values (m x k or k x m)
  n : desired number of values
  
  >sequence("x[n-1]+x[n-2]",[1,1],10)
  [1,  1,  2,  3,  5,  8,  13,  21,  34,  55]
  >sequence("n*x[n-1]",1,9)
  [1,  2,  6,  24,  120,  720,  5040,  40320,  362880]
  >A=[1,2;3,4]; sequence("A.x[,n-1]",[1;1],5)
  1          3         17         91        489
  1          7         37        199       1069
  >matrixpower(A,4).[1;1] // more efficiently
  489
  1069
  >function f(x,n) := mean(x[n-4:n-1])
  >plot2d(sequence("f",1:4,20)):
  
  See: 
iterate (Numerical Algorithms),
call (Euler Core)
function steffenson (f$:call, x:scalar, n:natural=0, eps=none)
  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 error 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 or a call collection of a variable x.
  
  f$ : function or expression of a scalar x
  x : start point
  n : optional number of iterations
  
  See: 
iterate (Numerical Algorithms),
call (Euler Core)
function nsteffenson (f$:call, x0:scalar, n:natural, eps=none)
  Does n steps of the Steffenson operator, starting from x0.
  
  Works like "steffenson", but returns the history of the iteration.
  
  See: 
niterate (Numerical Algorithms),
call (Euler Core)
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.
  
  >v=iterate("cos",1,40); longest v[-1]
  0.7390851699445544
  >longest aitkens(v)[-1]
  0.7390851332151597
  >longest iterate("cos",1)
  0.7390851332157367
  
  See: 
steffenson (Numerical Algorithms)

Newton Methods

function newton (f$:call, df$:call, x:complex scalar,
    y:number=0, eps=none)
  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 or call collection 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
  
  >newton("x^2","2*x",1,y=2)
  1.41421356237
  >function f(x) &= x^3-x+exp(x);
  >function df(x) &= diff(f(x),x);
  >newton(f,df,1), f(%)
  -1.13320353316
  0
  >newton(f,df,I), f(%)
  0.351537+0.802637i
  0+0i
  
  >function f(v) &= [v[1]+v[2]-1,v[1]^2+v[2]^2-10]
  >function Df(v) &= jacobian(f(v),[v[1],v[2]])
  >newton2("f","Df",[4,3]), f(%)
  
  See: 
secant (Numerical Algorithms),
bisect (Numerical Algorithms),
newton2 (Numerical Algorithms),
inewton (Numerical Algorithms),
mxminewton (Maxima Functions for Euler),
call (Euler Core)
function newton2 (f$:call, df$:call, x:numerical, eps=none)
  Newton method to solve a system of equations.
  
  The multidimensional Newton method can solve the equation f(v)=0
  for vectors v. Here, f must map row vectors to row vectors. This
  function needs the Jacobian matrix of f$, which is provided through
  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: 1 x n or n x 1, result 1 x n or n x 1)
  df$ : function df$(v) (v: 1 x n or n x 1, result: n x n)
  x : start point (1xn or nx1)
  
  >function f(x) := [x[1]+x[2]-1,x[1]^2+x[2]^2-10]
  >function Df(x) := [1,1;2*x[1],2*x[2]]
  >newton2("f","Df",[4,3]), f(%)
  [2.67945,  -1.67945]
  [0,  0]
  >function f([x,y]) &= [exp(x)-y,y-x-2];
  >function Df([x,y]) &= jacobian(f(x,y),[x,y]);
  >newton2("f","Df",[1,1]), f(%)
  [1.14619,  3.14619]
  [0,  0]
  
  See: 
broyden (Numerical Algorithms),
neldermin (Numerical Algorithms),
inewton2 (Numerical Algorithms),
call (Euler Core)
function broyden (f$:call, xstart:real, A:real=none, eps=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 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.
  
  f$ can be an expression or a call collection.
  
  To change the accuracy, you can specify an optional eps=...
  
  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),
call (Euler Core)
function nbroyden (f$:call, xstart:real vector, nmax:natural,
    A:real=none, eps=none)
  Do nmax steps of the Broyden algorithm.
  
  This function works like "broyden", but returns the steps of the
  algorithm. But it does at most nmax steps.
  
  f$ : function of one row vector, expression or call collection
  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.
  
  See: 
call (Euler Core)

Numerical Differentiation

DifMatrix:=%setupdif(5);
function diff (f$:call, 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 inherently somewhat inaccurate for
  general functions. To get a good approximation, the first
  derivative uses 4 evaluations of the function. There is a more
  accurate function diffc() for functions, which are analytic and
  real valued on the real line.
  
  f is either a function in a scalar, or an expression or call
  collection in x. Additional parameters after a semicolon are passed
  to f.
  
  f$ : function or expression of a scalar x
  x : the point, where the derivative should be computed
  n : degree of derivative
  h : optional delta
  
  >diff("x^x",2)
  6.77258872224
  >&diff(x^x,x)(2)
  6.77258872224
  
  See: 
diffc (Numerical Algorithms),
call (Euler Core)
function diffc (f$:call, x:numerical, h:number=epsilon())
  Computes the first derivative for real analytic functions
  
  This uses the imaginary part of a f(x+ih)/h. Thus it must be
  possible to evaluate f into the complex plane.
  
  f$ : function or expression. The function must evaluate for complex
  values.
  h : step size
  
  See: 
diff (Maxima Documentation),
call (Euler Core)

FFT, Folding, Filters

Most functions here are built-in functions. For an introduction to the Fourier Transform, see the following tutorial.

function comment overwrite fft (v:vector)
  Fast Fourier Transform of v
  
  v : Real or complex row vector. It should have a number of columns
  with many low prime factors, like 2^n. This is the same as
  evaluating a polynomial with coefficients v in all roots of unity
  simultaneously, but much faster.
  
  v : If v is a matrix, the function returns the two dimensional
  FFT. In this case, the number of rows and columns of v must be a
  power of 2.
  
  >sec=1.2; t=soundsec(sec); s=sin(440*t)+cos(660*t)/2;
  >i=1:2000; plot2d(i/sec,abs(fft(s))[i]):
  
  See: 
fht (Numerical Algorithms)
function comment overwrite ifft (v:vector)
  Inverse Fast Fourier Transform of v
  
  The inverse of fft.
  
  See: 
fft (Maxima Documentation)
function comment overwrite fht (v:real vector)
  Fast Hartley Transform of v
  
  The Hartley transform is similar to the Fourier Transform, but
  works from real vectors to real vectors. It can be used for sound
  analysis in the same way the Fourier transform can be used.
  
  >sec=1.2; t=soundsec(sec); s=sin(440*t)+cos(660*t)/2;
  >i=1:2000; plot2d(i/sec,fht(s)[i]):
  
  Algorithm from AlgLib.
  
  See: 
fft (Numerical Algorithms),
fft (Maxima Documentation)
function comment overwrite ifht (v:real vector)
  Inverse Fast Hartley Transform of v
  
  The inverse of fht().
  
  Algorithm from AlgLib.
  
function fftfold (v:vector, w:vector)
  fold v with w using the 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)
  
  >fold(1:10,[1,2,1]/4)
  [2,  3,  4,  5,  6,  7,  8,  9]
  >tail(fftfold(1:10,[1,2,1]/4),3)
  [2,  3,  4,  5,  6,  7,  8,  9]
  
function comment overwrite fold (v:complex, w:complex)
  fold v with w
  
  v and w can be real or complex vectors or matrices.
  
  In the result R, R[i,j] is the total sum of the pointwise
  multiplication of the submatrix of A with A[i,j] in the upper left
  corner and B. I.e.
  
  The result R is of size (cols(A)-cols(B)+1,rows(A)-rows(B)+1).
  Consequently B must be smaller in size than A.
  
  In the example, folding with [1/2,1/2] takes the average of two
  neighboring elements of 1:10.
  
  >fold(1:10,[0.5,0.5])
  [ 1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5 ]
  
  For a matrix, folding with [1/4,1/4;1/4,1/4] takes the average of
  2x2-sub-matrices.
  
  Folding with [1,-1] takes the difference of two elements, and is
  equivalent to differences.
  
  See: 
fftfold (Numerical Algorithms)
function overwrite filter (b:complex vector, a:complex vector, ..
    x:complex vector, y:complex vector=none, zeros:integer=true)
  Apply a filter with a and b to start values x and y
  
  Computes y[n] recursively using
  
  a[1]*y[n]+...+a[m]*y[n-m+1] = b[1]*x[n]+...+b[k]*x[n-k+1]
  
  a[1] must not be zero, of course.
  
  The start values on the right hand side are either 0 or the first k
  values of x, depending on the flag zeros (true by default). The
  start values on the left hand side are either 0 (the default), or
  the values of the parameter y.
  
  The size of a must exceed the size of x by 1.
  All input vectors must be real or complex row vectors or scalars.
  
  See: 
fold (Numerical Algorithms),
fftfold (Numerical Algorithms)

Interpolation

Aliases to function in Euler Core

function interpolate (x:vector, y:vector)
  Interpolate the values in y at x
  
  Alias to interp. Return the vector of divided differences.
  
  See: 
interp (Euler Core),
evaldivdif (Numerical Algorithms)
function evaldivdif (t:numerical, d:vector, x:vector)
  Evaluate the divided differences d at t
  
  Works like interpval with other order of parameters. The first
  parameter contains the points, where the polynomial should be
  evaluated.
  
  Alias to interpval with t-argument first.
  
  Example:
  >xp=0:0.1:1; yp=sqrt(xp);
  >dp=divdif(xp,yp);
  >plot2d("sqrt(x)-evaldivdif(x,dp,xp)",0,1):
  
  See: 
divdif (Euler Core)
function evalpoly (t:numerical, p:vector)
  evaluates a polynomial p in t.
  
  Polynomials are stored as a row vector with the constant
  coefficient in front. The function works for complex polynomials
  too.
  
  Alias to polyval with t-argument first.
  
  See: 
polyval (Euler Core)

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),
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
  simultaneous 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 approximation
  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).
  
  >x=equispace(-1,1,500); y=abs(x);
  >{t,d,h,r}=remez(x,y,10);
  >plot2d(x,interpval(t,d,x)-y):
  
  See: 
interpval (Euler Core),
polyfit (Linear Algebra)

Nonlinear Optimization

function fmin (f$:call, a:real scalar, b:real scalar=none,
    d=0.01, dmax=1e20, eps=none)
  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 evaluations for full accuracy.
  If no interval is given (b==none), the function searches left or
  right of a with initial step size d doubling d in each step until
  the function gets larger. Then the usual method is started.
  
  f$ is either an expression or call collection 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.
  
  >fmin("x^x",0,1), 1/E
  0.367879431154
  0.367879441171
  
  >function f(x) := x^x
  >f(fmin("f",0,1))
  0.692200627555
  >fmin({{"x^2+x-a",a=2}},-2,2)
  -0.500000015805
  
  >function f(x,a) := x^(a*x)
  >f(fmin("f",0,1;1.1),1.1)
  0.667198694057
  
  See: 
fmax (Numerical Algorithms),
call (Euler Core)
function fmax (f$:call, a:real scalar, b:real scalar=none,
    d=0.01, dmax=1e20, eps=none)
  Compute the maximum of the concave function f in [a,b].
  
  Works like fmin().
  
  See: 
fmin (Numerical Algorithms)
function fextrema (f$:call, a:real scalar, b:real scalar,
    n:integer=100, eps=none)
  Find all internal extremal points of f in [a,b].
  
  f may be an expression or call collection in x or a function.
  Additional parameters after a semicolon are passed to a function f.
  f must vectorize to its first argument.
  
  n : number of sub-intervals to be scanned
  eps : epsilon for the search of the maximum
  
  Returns {minima,maxima} (vectors, possibly of length 0)
  
  See: 
fmax (Numerical Algorithms),
fmin (Numerical Algorithms),
fzeros (Numerical Algorithms)
function fzeros (f$:call, a:real scalar, b:real scalar,
    n:integer=100, eps=none)
  Find all zeros of f in [a,b]
  
  f may be an expression or a call collection in x or a function.
  Additional parameters after a semicolon are passed to a function f.
  f must vectorize to its first argument
  
  The function uses fextrema() to find all local extreme values and
  uses bisectin() to find the zeros between sign changes.
  
  n : number of sub-intervals to be scanned eps : epsilon for
  the search of the maximum
  
  See: 
fextrema (Numerical Algorithms)
function brentmin (f$:call, a:real scalar, d=0.1, eps=epsilon())
  Compute the minimum of f around a with Brent's method.
  
  This function is no longer recommended. Use fmin() instead.
  
  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 or call collection in x, or a function in f.
  
  Return the point of minimum.
  
  >function f(x,a) := x^2+x-a
  >brentmin({{"f",2}},1;2)
  -0.499999992146
  
  See: 
fmin (Numerical Algorithms),
fmax (Numerical Algorithms),
call (Euler Core)
function neldermin (f$:call, 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 be applied to
  solve a system of equations by minimizing the norm of the errors.
  
  This function cannot be used for one dimension. Use fmin() instead
  with a proper interval to search for the minimum.
  
  f must be function of a row or column vector x, returning a real
  value. Additional parameters after the semicolon are passed to f.
  f can also be an expression or a call collection 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).
  
  >function f([x,y],a) &= x^2+y^2+exp(-x)+x*y
  
  2          - x    2
  y  + x y + E    + x
  
  >xmin=nelder({{"f",3}},[1,2])
  [0.432563,  -0.216282]
  >plot2d(&f(x,y,3),>contour); plot2d(xmin[1],xmin[2],>points,>add):
  
  See: 
fmin (Numerical Algorithms),
nelder (Euler Core),
call (Euler Core)
function overwrite nelder (f$:call, v:real, d=0.1, eps=epsilon())
  Compute the minimum of f around v with the Nelder-Mead method.
  
  For more explanations and examples, see neldermin().
  
  See: 
neldermin (Numerical Algorithms),
call (Euler Core)
function nlfit (f$:call, Df$:call, v:real)
  Minimizes f(x) from start vector v.
  
  This method is named after Levenberg-Marquardt. It minimizes a
  function from n variables to m variables (m>n) by minimizing the
  linear tangent function. If the norm gets larger it decreases the
  step size, until the step size gets closer than epsilon.
  
  While f can be an expression or a call collection, it is easier to
  use symbolic functions as in the example below.
  
  f(x) maps 1xn vectors to 1xm vectors (m>n)
  Df(x) is the Jacobian of f.
  
  >function f([x,y]) &= [x,y,x^2+2*y^2+1];
  >function Df([x,y]) &= jacobian(f(x,y),[x,y]);
  >nlfit("f","Df",[0.2,0.1])
  [-9.19403e-011,  -1.07327e-010]
  >function h(v) := norm(f(v))
  >nelder("h",[3,1]) // for comparison
  [6.33484e-007,  -1.9733e-007]
  
  See: 
nelder (Euler Core),
nelder (Numerical Algorithms),
call (Euler Core)
function modelfit (f$:call, param:vector, x:real, y:vector,
    powell=false, p=2, d=1e-2, eps=epsilon(), tries=2)
  Fit a model f$ to data x,y with initial guess p
  
  This uses the Powell or the Nelder-Mead algorithm to fit a
  non-linear function to data x[i],y[1].
  
  f$ : The model function f$(x,p).
  The function must vectorize to x. x will be a n x m vector for m
  data of n variables. The result of f$(x,p) must be a 1 x m vector.
  f$ can also be an expression or a call collection.
  
  x,y : The data to be used. x can be m x n, and y a 1 x n vector.
  param : Initial guess of the model parameter.
  x,y : Row vector of data.
  powell : If false, Nelder-Mead will be used.
  
  The following parameters are relevant for the Nelder-Mead method.
  Especially, decreasing d may help for difficult models.
  
  d : initial step size
  eps : accuracy
  tries : number of new starts
  
  >function model(x,p) := p[1]*cos(p[2]*x)+p[2]*sin(p[1]*x);
  >xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
  >ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143, ...
  >  1.91091,0.919576,-0.730975,-1.42001];
  >modelfit("model",[1,0.2],xdata,ydata)
  
  See: 
neldermin (Numerical Algorithms),
powellmin (Numerical Algorithms),
call (Euler Core)
function powellmin (f$:call, x:real vector)
  Minimize the function f with Powell's method
  
  f$ : function mapping a row vector to a real number
  x : start vector
  
  This function uses Powell's algorithm to find the minimum of a
  convex function. Make sure that the function is defined everywhere.
  You should return a large value outside the proper range of
  definition of f$. Alternatively, setting "errors off" will often
  work.
  
  Additional semicolon parameters are passed to f$. f$ can be a
  call collection or an expression in the vector xn.
  
  >function f([x,y]) := x^2+x*y+y^2+x
  >powellmin("f",[0,0])
  [-0.666667,  0.333333]
  >function f(x,a) := x^2-a*x
  >powellmin({{"f",3}},1)
  1.50000000147
  >function f(x) := log(x)*x
  >errors off; powellmin("f",1)
  0.367879439331
  >errors on;
  
  See: 
neldermin (Numerical Algorithms),
errors (Euler Core),
errors (Maxima Documentation)

Bezier Curves

function bezier (x,y)
  Compute the Bezier curve to points y in x
  
  The Bezier curve is the sum of y_i*B_i(t). The points y
  form the control polygon for the curve. The curve interpolates
  only the first and last point of y.
  
  In the following example, we take a square as the control
  polygon. Note that t runs from 0 to 1 along the curve.
  
  >y=[0,0,1,1;0,1,1,0];
  >t=linspace(0,1,500);
  >plot2d(y[1],y[2],>addpoints);
  >p=bezier(t,y);
  >plot2d(p[1],p[2],>add):
  
  See: 
spline (Numerical Algorithms)

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 map splineval (t:number; x:vector, y:vector, h:vector)
  Evaluates the cubic interpolation spline for [x(i),y(i)]
  
  The function needs the second derivatives h[i] at t[i].
  
function map 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 map 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, depending on the value of the parameter "constant".
  
  >xp=1:10; yp=intrandom(1,10,3);
  >plot2d("linsplineval(x,xp,yp,<constant)",0,11);
  >plot2d(xp,yp,>points,>add);
function map 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)
function nspline (x:real, k:column integer, n:scalar integer,
    t:real vector)
  Compute the normalized B-spline on t[k,k+n+1] in x.
  
  This spline is the normalized B-spline on the interval
  (t[k],t[k+n+1]) and has knots t[k] to t[k+n+1]. The polynomial
  degree is n. The spline is evaluated in the point x.
  
  This function maps to vectors x, but not to k and n.
  
  Examples:
  
  >plot2d("nspline";1,2,[1,2,3,4],a=1,b=4):
  >plot2d("nspline";(1:3)',2,[1,2,3,4,5,6],a=1,b=6):
  

Adaptive Evaluation

There are many adaptive functions in EMT for integration and for differential equations. Moreover, the plot2d() function evaluates a function adaptively by default using adaptiveeval() and adaptiveevalone(). The function adaptive() is for adaptive evaluation of a general function.

function adaptive (f$:call, a:number, b:number,
    eps=0.01, dmin=1e-5, dmax=1e5, dinit=0.1, factor=2)
  Compute f$(x) for x in [a,b] with adaptive step size.
  
  The adaptive evaluation takes step sizes of different sizes to be
  able to return the values of the function with spacing not
  exceeding a give accuracy. The function is evaluated on an
  interval, with real or complex scalar or vector values.
  
  f$ is an expressions in x, a function or a call collection. The
  function must return a real or complex scalar or a column vector.
  
  eps is the target accuracy.
  a,b are the interval limits
  dmin and dmax are the minimal and maximal step size.
  dinit is the initial step size.
  factor is the factor the step size is changed with.
  
  Returns {x,y}
  x are the points of evaluation, and y is a row vector of values (or
  a matrix with columns of vector values).
  
  >{x,y}=adaptive("sqrt(x)*exp(-x)",0,10); plot2d(x,y):
  >plot2d(differences(x),>logplot):
  
  >{x,y}=adaptive("x*sin(1/x)",1,epsilon); plot2d(x,y):
  
  >function f(t:number,a:vector) := sum(exp(a*I*t));
  >{t,z}=adaptive({{"f",[1,2,5]}},0,2*pi); plot2d(z,r=3):
  
  See: 
adaptiveeval (Plot Functions)

Numerical Integration

Numerical integration in EMT can be done with various methods. The most accurate methods are the adaptive methods adaptivegauss() or alintegrate(). Less accurate, but fast methods are the Gauss quadrature or the LSODA integration. The function integrate defaults to the adaptive Gauss method.

function args simpson (f$:call, a:number, b:number,
    n:natural=50, maps:integer=0, allvalues=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.
  Additional parameters 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 (Numerical Algorithms)
function simpsonfactors (n)
  Returns the vector 1,4,2,...,2,4,1.
  
function romberg (f$:call, a:number, b:number, m:natural=10,
    maps:integer=0, eps=none)
  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.
  Additional parameters 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 expression needs to be mapped. Functions are always
  mapped.
  
  >romberg("x^x",0,1)
  0.783430510714
  
  See: 
simpson (Numerical Algorithms),
gauss (Numerical Algorithms),
integrate (Numerical Algorithms),
integrate (Maxima Documentation)
function adaptiveint (f$:call, a:real scalar, b:real scalar, ..
    eps=epsilon(), steps=10)
  Returns the integral from a to b with the adaptive Runge method.
  
  f is an expression in x, which must not contain y, or a function
  in x. Make sure, that f is computed at each point with relative
  error eps or less. Additional parameters after a semicolon are
  passed to a function f.
  
  f : function or expression in x
  a,b : interval bounds
  eps : accuracy of the adaptive process
  steps : hint for a good step size
  
  >adaptiveint("x^x",0,1)
  0.783430510714
  
  See: 
adaptiverunge (Numerical Algorithms),
integrate (Numerical Algorithms),
integrate (Maxima Documentation)
function adaptiveintlsoda (f$:call, a:real scalar, b:real scalar, ..
    eps=epsilon(), steps=10)
  Returns the integral from a to b with the adaptive Runge method.
  
  f is an expression in x, which must not contain y, or a function
  in x. Make sure, that f is computed at each point with relative
  error eps or less. Additional parameters after a semicolon are
  passed to a function f.
  
  f : function or expression in x
  a,b : interval bounds
  eps : accuracy of the adaptive process
  steps : hint for a good step size
  
  See: 
adaptiverunge (Numerical Algorithms),
integrate (Numerical Algorithms),
integrate (Maxima Documentation)
function overwrite alintegrate (f$:call, a:number, b:number)
  Adaptive integration
  
  Algorithm  from AlgLib similar to adaptivegauss().
  
  See: 
integrate (Maxima Documentation)
function overwrite alsingular (f$:call, a:number, b:number,
    alpha:number, beta:number)
  Adaptive integration
  
  Algorithm from AlgLib to integrate a function with singularities
  at the boundaries. The singularities must be of type (x-a)^alpha
  and (b-x)^beta. Either alpha or beta can be 0, which denotes a
  smooth function.
  
  >alsingular("1/sqrt(x)",0,1,-0.5,0)
  2
  
  See: 
integrate (Numerical Algorithms),
integrate (Maxima Documentation),
alintegrate (Numerical Algorithms)
function integrate (f$:call, a:real scalar, b:real scalar, ..
    eps=epsilon(), steps=10, method:integer=0, fast:integer=0,
    maps:integer=1)
  Integrates f from a to b with the adaptive Runge method.
  
  Calls an integration method, depending on the method parameter. The
  default is the Gauss integration with 10 subintervals. Available is
  the adaptive Runge, Gauss and the LSODA method. Note, that the
  fastest algorithm is the Gauss algorithm, which is the method of
  choice.
  
  method :
  0=adaptive Gauss, (needs vectorized function)
  1=adaptive Runge,
  2=lsoda,
  3=AlgLib
  
  fast : Simple Gauss with one sub-interval (needs vectorized
  function)
  
  >integrate("x^x",0,1)
  0.783430510712
  
  See: 
mxmiint (Maxima Functions for Euler),
gauss (Numerical Algorithms),
adaptivegauss (Numerical Algorithms),
alintegrate (Numerical Algorithms),
lsoda (Numerical Algorithms)

Gauss-Quadrature

function legendre (n:index)
  Compute the coefficients of the n-th Legendre polynomial.
  
  This is used to compute the Gauss coefficients.
  
  See: 
gauss (Numerical Algorithms)
function gaussparam (n:index)
  Returns the knots and alphas of gauss integration at n points in [-1,1].
  
  Returns {gaussz,gaussa}.
  
function gauss (f$:call, a:real scalar, b:real scalar, n:index=1, ..
    maps:integer=1, xn=gaussz, al=gaussa)
  Gauss integration with 10 knots and n subintervals.
  
  This function is exact for polynomials up to degree 19, even for
  one sub-interval. For other functions, it may be necessary to
  specify n subintervals.
  
  f is an expression in x, or a function in x. Additional parameters
  after the semicolon are passed to f.
  
  maps : Vectorize an expression to the arguments. On by default.
  Turn off for a bit more performance.
  
  See: 
gauss5 (Numerical Algorithms),
simpson (Numerical Algorithms),
integrate (Numerical Algorithms),
integrate (Maxima Documentation),
romberg (Numerical Algorithms),
romberg (Maxima Documentation)
function igauss (expr$:call, a:real, b:real, n:index=1, ..
    dexpr20$:string=none, dexpr20est:real=none, xn=gaussz, al=gaussa)
  Compute an interval inclusion for an integral
  
  This works only for expressions, unless you provide an estimate for
  the 20-th derivative. The function will compute this derivative
  with Maxima otherwise.
  
  Optional:
  
  dexpr20: Expression or function for the 20-th derivative
  dexpr20est: Interval inclusion of the 20-th derivative
  
  >igauss("exp(-x^2)",0,1)
  ~0.7468241328124256,0.7468241328124285~
  
  See: 
gauss5 (Numerical Algorithms)
function gauss5 (f$:call, a:real scalar, b:real scalar, n:index=1, ..
    maps:integer=0, xn=gauss5z, al=gauss5a)
  Gauss integration with 4 knots and n subintervals
  
  This is exact for polynomials to degree 7 (even for n=1).
  The results are usually much better than the Simpson method
  with only twice as many function evaluations.
  
  See: 
simpson (Numerical Algorithms),
integrate (Numerical Algorithms),
integrate (Maxima Documentation),
romberg (Numerical Algorithms),
romberg (Maxima Documentation)
function gaussfxy (f$:call, a:real scalar, b:real scalar, ..
    c:real scalar, d:real scalar, n:index=1)
  Computes the double integral of f on [a,b]x[c,d].
  
  This function uses Gauss integration for the inner and the outer
  integral. The function is exact for polynomials in x and y up to
  degree 19.
  
  f$ : a function f(x,y) or an expression of x and y.
  
  >gaussfxy("exp(x)*y+y^2*x",0,1,0,2)
  4.76989699025
  >&integrate(integrate(exp(x)*y+y^2*x,x,0,1),y,0,2)() // exact!
  4.76989699025
  
function adaptivegauss (f$:call, a:real scalar, b:real scalar, ..
    eps=epsilon*100, n=1, maps:integer=1)
  Adaptive Gauss integral of f on [a,b].
  
  This function divides the interval into subintervals, and tries a
  Gauss integration with 1 and 2 intervals on each sub-interval. If
  both do not agree good enough the step size is divided by 2. the
  algorithm tries to double the step size in each step..
  
  n : default number of subintervals
  eps : default accuracy
  maps : Vectorizes the call to f$. Turn off for more performance.
  
  Examples:
  adaptivegauss("sqrt(x)",0,1)
  

Differential Equations

The LSODA method should be the method of choice. It is a very efficient adaptive method for ordinary differential equations, which takes also care of stiff cases. But the Heun and Runge methods are available too, including an adaptive Runge method, which can be used for adaptive integration.

Vector fields can be plotted using the vectorfield() function in plot2d.

There are also the interval algorithms mxmidlg() and mxmiint().

Numerical Solutions of Differential Equations

function heun (f$:string, t:real vector, y0:real vector)
  Solves the differential equation y'=f(x,y) at x=t with y(t[1])=y0.
  
  This function works also for systems of n differential equations.
  
  f must be a function with two parameters f(x,y), where x is scalar,
  and y is scalar or a 1xn row vector in the case of a system of
  differential equations. f must return a scalar, or a 1xn row vector
  respectively. Additional parameters for "heun" after a semicolon
  are passed to a function f. Alternatively, f can be an expression
  in x and y.
  
  y0 is the initial value y(t[1]), a scalar, or a 1xn row vector.
  t is the 1xm row vector of points, where the differential equation
  will be solved.
  
  f : function or expression in x and y
  t : row vector of points, where y(t) should be computed
  y0 : initial value, scalar or row vector
  
  The function returns a 1xm row vector of y, or a nxm matrix in the
  case of a system of n equations.
  
  >x=0:0.01:2; plot2d(x,heun("-2*x*y",x,1)):
  >function f(x,y,a) &= -a*x*y;
  >heun("f",x,1;2)[-1]
  0.0183155649842
  >exp(-4)
  0.0183156388887
  
  See: 
runge (Numerical Algorithms)
function runge (f$:string, t:real vector, y0:real, ..
    steps:index=1)
  Solves the differential equation y'=f(x,y) at x=t with y(t[1])=y0.
  
  This function works also for systems of n differential equations.
  
  f must be a function with two parameters f(x,y), where x is scalar,
  and y is scalar or a vector in the case of a system of differential
  equations. f must return a scalar, or a vector respectively.
  Additional parameters for "runge" after a semicolon are passed to a
  function f. Alternatively, f can be an expression in x and y.
  
  y0 is the initial value y(t[1]), a scalar, or a 1xn row vector.
  t is the 1xm vector of points, where the differential equation will
  be solved.
  
  steps are optional intermediate steps between the t[i]. This
  parameter allows to take extra steps, which will not be stored into
  the solution.
  
  Note that f can work with row or column vectors, but the initial
  value must be of this form too. The return matrix will always
  consist of column vectors.
  
  f : function or expression in x and y
  t : row vector of points, where y(t) should be computed
  y0 : initial value, scalar or row vector
  steps : additional steps between the points of x
  
  Returns a row vector y(t[i]) for the elements of t. Im case of a
  system, the function returns a matrix with columns y(t[i]).
  
  >x=0:0.01:2; plot2d(x,runge("-2*x*y",x,1)):
  >function f(x,y,a) &= -a*x*y;
  >runge("f",x,1;2)[-1]
  0.018315639424
  >exp(-4)
  0.0183156388887
  
  See: 
heun (Numerical Algorithms),
adaptiverunge (Numerical Algorithms)
function adaptiverunge (f$:call, x:real vector, y0:real, ..
    eps=epsilon(), initialstep=0.1)
  Solves the differential equation y'=f(x,y) at x=t with y(t[1])=y0
  adaptively.
  
  The function uses an adaptive step size between x[i] and x[i+1] to
  guarantee an accuracy eps. By default the accuracy is the default
  epsilon, but eps=... can be specified as an extra parameter.
  
  This function works also for systems of n differential equations.
  
  f must be a function with two parameters f(x,y), where x is scalar,
  and y is scalar or a vector in the case of a system of differential
  equations. f must return a scalar, or a vector respectively.
  Additional parameters for "adaptiverunge" after a semicolon are
  passed to a function f. Alternatively, f can be an expression in x
  and y.
  
  y0 is the initial value y(t[1]), a scalar or a vector. t is the 1xm
  vector of points, where the differential equation will be solved.
  
  Note that f can work with row or column vectors, but the initial
  value must be of this form too. The return matrix will always
  consist of column vectors.
  
  f : function or expression in x and y
  t : row vector of points, where y(t) should be computed
  y0 : initial value, scalar or row vector
  eps : accuracy of the adaptive method
  initialstep : initial step size
  
  Returns a row vector y(t[i]) for the elements of t. Im case of a
  system, the function returns a matrix with columns y(t[i]).
  
  >longest adaptiverunge("-2*x*y",[0,2],1)
  1     0.01831563888880886
  >longest E^-4
  0.01831563888873419
  
  See: 
heun (Numerical Algorithms)
function ode (f$:call, t:real vector, y:real, ..
    eps:real=epsilon(), warning:integer=false, reset:integer=false)
  Solves the differential equation y'=f(x,y) at x=t with y(t[1])=y0.
  
  This function works also for systems of n differential equations.
  The algorithm used is the LSODA algorithm for stiff equations.
  This algorithm switches between stiff and normal case
  automatically.
  
  The LSODA algorithm is based on work by Linda R. Petzold  and
  Alan C. Hindmarsh, Livermore National Laboratory. The EMT version
  is based on a C source by Heng Li, MIT.
  
  f must be a function with two parameters f(x,y), where x is scalar,
  and y is scalar or a vector in the case of a system of differential
  equations. f must return a scalar, or a vector respectively.
  Additional parameters for "lsoda" after a semicolon are passed to a
  function f. Alternatively, f can be an expression or a call
  collection in x and y.
  
  This is an adaptive algorithm. But for some functions, it is
  necessary to add intermediate points, especially if the function is
  close to 0 over wide subintervals. Otherwise, it has been observed
  that the algorithm fails. If you are sure that the problem does not
  fall under the problematic category, you can leave reset=false, which
  will interpolate intermediate points instead of resetting the
  algorithm every time. This is a huge benefit in terms of function
  evaluations. By default, reset is false.
  
  y0 is the initial value y(t[1]), a scalar, or a 1xn row vector.
  t is the 1xm vector of points, where the differential equation will
  be solved.
  
  Note that f can work with row or column vectors, but the initial
  value must be of this form too. The return matrix will always
  consist of column vectors.
  
  f : function or expression in x and y
  t : row vector of points, where y(t) should be computed
  y0 : initial value, scalar or row vector
  warning : toggles warnings from lsoda (but not errors)
  reset : reset the integration between iterations
  
  Returns a row vector y(t[i]) for the elements of t. Im case of a
  system, the function returns a matrix with columns y(t[i]).
  
  >function f(x,y,a) &= -a*x*y;
  >ode({{"f",2}},x,1)[-1]
  0.0183156388888
  >exp(-4)
  0.0183156388887
  
  >ode("x*y",[0:0.1:1],1)[-1], exp(1/2) // solve y'=x*y
  1.64872127078
  1.6487212707
  
  >function f(x,y) := -x*y^2/(y^2+1); // y'=-xy/(y^2+1)
  >t=linspace(0,5,1000); plot2d(t,ode("f",t,1)):
  
  >sol &= ode2('diff(y,x)=-x*y^2/(1+y^2),y,x) // check with Maxima
  >&solve(sol with %c=0,y)[2], plot2d(&rhs(%),0,5):
  
  >function f(x,y) := [y[2],-10*sin(y[1])]; // y''=-10*sin(y)
  >t=linspace(0,2pi,1000); plot2d(t,ode("f",t,[0,1])[1]):
  
  See: 
heun (Numerical Algorithms),
runge (Numerical Algorithms),
adaptiverunge (Numerical Algorithms),
lsoda (Numerical Algorithms),
call (Euler Core)
function overwrite lsoda (f$:call, t:real vector, y:real, ..
    eps:real=epsilon(), warning:integer=true, reset:integer=false)
  Solves the differential equation y'=f(x,y) at x=t with y(t[1])=y0.
  
  See: 
ode (Numerical Algorithms)

Sparse Matrices

Euler has support for a compressed format for thin matrices. This format stores non-zero elements only. Operations on such matrices are typically very much faster. For examples see the following introduction.

To solve large, sparse equations the CG-method implemented in cgX() or cpxfit() is the method of choice. The function cpxfit() uses the normal equation to fit Ax-b. The Gauss-Seidel method seidelX() is an alternative for diagonal dominant matrices.

function comment cpx (A)
  Compress the matrix A
  
  The compressed matrix stores only the non-zero values of A. It
  consists of lines of the form (i,j,x), which means A[i,j]=x. Note
  that compressed matrices are a separate data type in Euler.
  
  >cpx(id(3))
  Compressed 3x3 matrix
  1     1                    1
  2     2                    1
  3     3                    1
  
  See: 
decpx (Numerical Algorithms)
function comment decpx (X)
  Decompress the compressed matrix X
  
  >C=cpxzeros(3,3)
  Compressed 3x3 matrix
  >decpx(cpxset(C,[1,1,0.5;2,3,0.7]))
  0.5             0             0
  0             0           0.7
  0             0             0
  
function comment cpxzeros ([n,m])
  Empty compressed matrix of size nxm
function comment cpxset (X,K)
  Set elements in the compressed matrix X to
  
  K is a matrix containing lines of type (i,j,x). Then X[i,j] will be
  set to x. If the element already exists in X, it is replaced by the
  new value.
  
  >function CI(n) := cpxset(cpxzeros(n,n),(1:n)'|(1:n)'|1)
  >decpx(CI(3))
  1             0             0
  0             1             0
  0             0             1
  
  See: 
decpx (Numerical Algorithms),
cpxget (Numerical Algorithms)
function comment cpxget (X)
  Get the elements of a compressed matrix X
  
  The return value is a matrix containing lines of type (i,j,x).
  
  >H=[1,1,0.1;2,3,0.7]
  1             1           0.1
  2             3           0.7
  >cpxget(cpxset(cpxzeros(3,3),H))
  1             1           0.1
  2             3           0.7
  
  See: 
cpxset (Numerical Algorithms)
function comment cpxmult (A,B)
  Multiply two compressed matrices
  
  This is a fast algorithm to multiply two compressed, sparse
  matrices. The result is a compressed matrix.
  

Gauss-Seidel for Sparse Systems

function comment cpxseidel (C,b,x,om)
  One step of the Gauss-Seidel algorithm for Cx=b.
  
  C must be a compressed matrix, x is the step to go from, and om is
  the relaxation parameter.
  
  See: 
seidelX (Numerical Algorithms),
cgX (Numerical Algorithms)
function seidelX (H:cpx, b:real column, x:column=0, ..
    om:real positive scalar=1.2, eps=none)
  Solve Hx=b using the Gauss-Seidel method for compressed matrices H.
  
  The Gauss-Seidel method  with Relaxation is an iterative method,
  which converges for all positive definite matrices. For large
  matrices, it may work well. However, the conjugate gradient method
  "cgX" is the method of choice.
  
  H must be diagonal dominant, at least not have 0 on the diagonal.
  om is the relaxation parameter between 1 and 2. x is start value
  (automatic if 0). It is possible to specify the accuracy with
  eps=...
  
  H : compressed matrix (nxm)
  b : column vector (mx1)
  om : optional relaxation coefficient
  
  Returns the solution x, a column vector.
  
  See: 
cgX (Numerical Algorithms),
seidel (Linear Algebra)
function cgX (H:cpx, b:real column, x0:real column=none,
    f:index=10, eps=none)
  Conjugate gradient method to solve Hx=b for compressed H.
  
  This is the method of choice for large, sparse matrices. In most
  cases, it will work well, fast, and accurate.
  
  H must be positive definite. Use cpxfit(), if it is not.
  
  The accuracy can be controlled with an additional parameter
  eps. The algorithm stops, when the error gets smaller then eps, or
  after f*n iterations, if the error gets larger. x0 is an optional
  start vector.
  
  H : compressed matrix (nxm)
  b : column vector (mx1)
  x0 : optional start point (mx1)
  f : number of steps, when the method should be restarted
  
  >X=cpxsetdiag(cpxzeros(1000,1000),0,2);
  >X=cpxsetdiag(cpxsetdiag(X,-1,1),1,1);
  >b=random(1000,1);
  >x=cgX(X,b);
  >totalmax(abs(cpxmult(X,x)-b))
  0
  
  See: 
cpxfit (Numerical Algorithms),
cg (Linear Algebra),
cgXnormal (Numerical Algorithms)
function cgXnormal (H:cpx, Ht:cpx, b:real column, x0:real column=none, ..
    f:index=10, eps=none)
  Conjugate gradient method to solve Ht.H.x=Ht.b for compressed H.
  
  This algorithm is used by cpxfit() to solve the normal equation
  H'Hx=H'b, which minimizes |Hx-b|, i.e., to find an optimal solution
  of an linear system with more equations than unknowns.
  
  Stops, when the error gets smaller then eps, or after f*n
  iterations, when the error gets larger. x0 is an optional start
  vector.
  
  H : compressed matrix (nxm)
  Ht : compressed matrix (mxn)
  b : column vector (mx1)
  x0 : optional start (mx1)
  f : number of steps, when the method should be restarted
  
  See: 
cpxfit (Numerical Algorithms)
function cpxfit (H:cpx, b:real column, f:index=10, eps=none)
  Minimize |Hx-b| for a compressed matrix H.
  
  This function uses the conjugate gradient method to solve the
  normal equation H'Hx=H'b for sparse compressed matrices H.
  
  H : compressed matrix (nxm)
  b : column vector (mx1)
  f : number of steps, when the method should be restarted
  
  >X=cpxsetdiag(cpxzeros(1000,1000),0,2);
  >X=cpxsetdiag(X,1,1);
  >b=random(1000,1);
  >x=cpxfit(X,b);
  >totalmax(abs(cpxmult(X,x)-b))
  0
  
  See: 
fit (Linear Algebra),
fitnormal (Linear Algebra),
svdsolve (Linear Algebra)
function cpxsetdiag (R:cpx, k:integer scalar, d:real vector)
  Set the k-th diagonal of the compressed matrix R to the value d.
  
  k=0 is the main diagonal, k=-1 the diagonal below, and k=1 the
  diagonal above.
  
  Note that this function does not change R, but returns a new matrix
  with the changes.
  
  >function CI(n) := cpxsetdiag(cpxzeros(n,n),0,1);
  >decpx(CI(3))
  1             0             0
  0             1             0
  0             0             1
  
  See: 
setdiag (Euler Core)
function cpxmultrows (c:real column, A:cpx)
  Multiply the i-th row of the compressed matrix A by c[i].
  
  Note that this function does not change R, but returns a new matrix
  with the changes.
  
  c : column vector (nx1)
  A : compressed matrix (nxm)
  
  >function CI(n) := cpxsetdiag(cpxzeros(n,n),0,1);
  >decpx(cpxmultrows((1:3)',CI(3)))
  1             0             0
  0             2             0
  0             0             3
  
  See: 
cpxsetrow (Numerical Algorithms)
function cpxsetrow (A:cpx, i:index, r:real vector)
  Set the i-th row of the compressed matrix A to r.
  
  Note that this function does not change R, but returns a new matrix
  with the changes.
  
  A : compressed matrix (nxm)
  i : index
  r : row vector (1xm)
  
function cpxsetcolumn (A:cpx, j:index, c:real column)
  Set the j-th column of the compressed matrix A to c.
  
  Note that this function does not change R, but returns a new matrix
  with the changes.
  
  A : compressed matrix (nxm)
  j : integer
  c : column vector (mx1)
  
  See: 
cpxsetrow (Numerical Algorithms)
function cpxsum (A:cpx)
  The sums of all rows of the compressed matrix A.
  
  See: 
sum (Euler Core),
sum (Maxima Documentation)

Incidence Matrices

function rectangleIncidenceX (n:index, m:index)
  Incidence matrix of a rectangle grid in compact form.
  
  The incidence matrix of a graph is the matrix H, such that H(i,j)
  contains 1, if node i is connected to node j. In this function, the
  graph consists of the points of a rectangle, and the edges connect
  adjacent points. The points in the rectangle are numbered row by
  row.
  
  Returns a compressed nm x nm matrix.
  
  See: 
cpxset (Numerical Algorithms)

Exact Computation

Euler has a long accumulator, which can compute the scalar product and the residuum of a linear equation exactly. On this basis, Euler implements a residuum iteration to solve linear systems. It can be used to invert matrices more exactly, or to evaluate polynomials.

Residuum iterations uses the long accumulator of Euler to compute an exact scalar product.

function xlgs (A:complex, b:complex, n:integer=20, eps=none)
  Compute a more exact solution of Ax=b using residuum iteration.
  
  You can specify the relative accuracy with eps=... This epsilon is
  used to determine, if the algorithm should stop.
  
  n : the maximal number of residual iterations.
  
  >H=hilbert(10);
  >longest totalmax(abs(xlgs(H,sum(H))-1))
  0
  >longest totalmax(abs(H\sum(H)-1))
  8.198655030211555e-005
  
  See: 
ilgs (Numerical Algorithms)
function xinv (A:complex, n:integer=20, eps=none)
  Compute the inverse of A using residuum iteration.
  
  You can specify the relative accuracy with eps=... as last
  parameter. Additionally, a maximal number of iteration n can be
  set.
  
  See: 
xlgs (Numerical Algorithms)
function xlusolve (A:complex, b:complex, n:integer=20, eps=none)
  Compute the solution of Ax=b for L- or U-matrices A.
  
  Works for lower triangle matrices with diagonal 1, or for upper
  triangle matrices. The function is just a faster version of xlgs.
  
  You can specify the relative accuracy with eps=... as last
  parameter. Additionally, a maximal number of iteration n can be
  set.
  
  See: 
xlgs (Numerical Algorithms),
lu (Linear Algebra),
lusolve (Linear Algebra)
function xpolyval (p:complex vector, t:complex, n:integer=20, eps=none)
  Evaluate the polynomial at values t using residuum iteration.
  
  Alias to xevalpoly().
  
  See: 
xevalpoly (Numerical Algorithms)
function xevalpoly (t:complex, p:complex vector,
    n:integer=20, eps=none)
  Evaluate the polynomial at values t using residuum iteration.
  
  You can specify the relative accuracy with eps=... as last
  parameter. Additionally, a maximal number of iteration n can be
  set.
  
  >p:=[-945804881,1753426039,-1083557822,223200658]; ...
  >t:=linspace(1.61801916,1.61801917,100); ...
  >plot2d(t-1.61801916,evalpoly(t,p)):
  >plot2d(t-1.61801916,xevalpoly(t,p,eps=1e-17)):
  
  See: 
xlgs (Numerical Algorithms)
function xeigenvalue (a,l,x=none)
  Returns an improved eigenvalue of A, starting with the approximation l.
  
  l must be close to a simple eigenvalue, and x close to an
  eigenvector, if x is not 0. This is the inverse iteration due to
  Wielandt.
  
  Returns the eigenvalue and the eigenvector.
  
  >H=hilbert(10);
  >lambda=jacobi(H)[1]
  255023613680
  >lambda=lambda+1000;
  >x=eigenspace(H,lambda);
  >norm(H.x-lambda*x)
  7445.4368515
  >{lambda,x}=xeigenvalue(H,lambda,x);
  >norm(H.x-lambda*x)
  3.7570410926e-005
  
  See: 
eigenvalues (Linear Algebra),
eigenvalues (Maxima Documentation)

Interval Solvers and Guaranteed Solutions

function ieval (f$:string, x:interval scalar, n:integer=10)
  Better evaluation of the expression in f for the interval x.
  
  Evaluating an expression directly assumes that all values in all
  intervals are combined. If the expression is a function of x, this
  is not the way we want to to evaluate the expression. Then we want
  to take the same point whenever x occurs in the expression. So
  the interval is split into subintervals for more accuracy, and then
  the usual evaluation takes place in each subinterval. The results
  are combined.
  
  n : number of subintervals.
  
  >expr &= x^3-x+1/x^2; plot2d(expr,0.8,0.9):
  
  >x=~0.8,0.9~; expr()
  ~0.84,1.5~
  >ieval(expr,~0.8,0.9~) // better
  ~0.99,1.3~
  >ieval(&diff(expr,x),~0.8,0.9~)<0 // decreasing function
  1
  >longest expr(~0.8~)||expr(~0.9~) // correct
  ~1.063,1.275~
  
  >function f(x) := @expr
  >ieval("f",x) // the same with functions
  ~0.99,1.3~
  
  See: 
mxmieval (Maxima Functions for Euler),
ievalder (Numerical Algorithms)
function ievalder (f$:string, fd$:string, xi:interval scalar, n:integer=10)
  Better evaluation of the expression in f for the interval x.
  
  The derivative is used to improve the interval accuracy. The
  interval is split into sub-intervals for more accuracy.
  
  See: 
mxmieval (Maxima Functions for Euler)
function idgl (f$:string, x:real vector, y0:interval scalar)
  Guaranteed inclusion of y'=f(t,y0;...) at points t with y(t[1])=y0.
  
  This is a quick inclusion for a differential equation, which avoids
  the use of any Taylor series. The inclusion is very coarse,
  however. The function uses a simple Euler one step method.
  
  The result is an interval vector of values.
  
  See: 
mxmidgl (Maxima Functions for Euler),
idglder (Numerical Algorithms)
function idglder (f$:string, fx$:string, fy$:string, x:real vector, ..
    y0:interval scalar)
  Inclusion for y'=f(t,y0;...) at t with y(t[1])=y0.
  
  This function needs the partial derivatives of f to x and y.
  
  The result is an interval vector of values. The result is a very
  coarse inclusion. For better results, see mxmidgl().
  
  f, fx and fy are functions in f(x,y), or expressions of x and y.
  Additional arguments are passed to the functions.
  
  >function f(x,y) &= -x*sin(y);
  >function fx(x,y) &= diff(f(x,y),x);
  >function fy(x,y) &= diff(f(x,y),y);
  >x=0:0.01:2;
  >idglder("f","fx","fy",x,~1~)[-1]
  ~0.137,0.158~
  >mxmidgl(&f(x,y),x,~1~)[-1]
  ~0.14759945743769,0.14759945743821~
  
  See: 
mxmidgl (Maxima Functions for Euler)
function isimpson (f$:string, der$:string, a:number, b:number, ..
    n:index=50)
  Interval Simpson integral of f from a to b.
  
  This function uses the Simpson method and its error estimate to get
  guaranteed inclusions of integrals. Other interval methods like
  igauss or mxmiint provide better estimates.
  
  f : expression (must map to arguments and work for intervals)
  der : expression for fourth derivative (like f)
  a,b : interval bounds
  n : number of subintervals
  
  >isimpson("exp(-x^2)",&diff(exp(-x^2),x,4),0,2)
  ~0.8820813879,0.8820813936~
  >igauss("exp(-x^2)",0,2)
  ~0.88208139037,0.88208139115~
  >mxmiint("exp(-x^2)",0,2)
  ~0.8820813907619,0.882081390763~
  
  See: 
mxmisimpson (Maxima Functions for Euler)
function ilgs (A:interval, b:interval, R="", steps=100)
  Guaranteed interval inclusion for the solution of A.x=b.
  
  This function uses an interval algorithm, and an exact residuum
  calculation. If the algorithm succeeds, the result is a guaranteed
  inclusion for the solution of the linear system. Note that the
  algorithm can only work for regular A, or interval matrices not
  containing singular A.
  
  A and b may be of interval or real type.
  
  The optional R should be a pseudo inverse to A.
  
  >H=hilbert(10); b=sum(H);
  >longest max(diameter(ilgs(H,b)'))
  1.695865670114927e-012
  >longest max(abs(xlgs(H,b)'-1))
  0
  >longest max(abs((H\b)'-1))
  8.198655030211555e-005
  
  See: 
xlgs (Numerical Algorithms)
function iinv (A:interval)
  Guaranteed interval inverse of the matrix A.
  
  See: 
inv (Linear Algebra),
xinv (Numerical Algorithms)
function ievalpoly (t:interval, p:interval vector)
  Guaranteed evaluation of a polynomial p(t).
  
  p contains the coefficients of a polynomial. Euler stores
  polynomials starting with the constant coefficient.
  
  See: 
polyval (Euler Core),
xpolyval (Numerical Algorithms)
function ipolyval (p:interval vector, t:interval)
  Guaranteed evaluation of a polynomial p(t).
  
  See: 
ievalpoly (Numerical Algorithms)
function ibisect (f:string, a:scalar, b:scalar=none, y:scalar=0)
  Interval bisection algorithm to solve f(x)=y
  
  >function f(x) &= x^2+x;
  >x=ibisect("f(x)",0,4,y=2)
  ~0.99999999999999967,1.0000000000000004~~
  >x=ibisect("f(x)",0,4,y=~1.99,2.01~)
  ~0.9958,1.005~~
  >function f(x,a) &= x^a-a^x;
  >ibisect("f",1,3;~1.999,2.001~) // a as semicolon parameter, y=0
  ~1.98,2.02~~
  
  See: 
inewton (Numerical Algorithms)
function ibis (f,a,b=none,y=0) ...
    if b==none then b=right(a); a=left(a); endif;


function inewton (f$:string, df$:string , xi: interval, yi:real scalar="", y=0)
  Guaranteed interval inclusion of the zero of f.
  
  df must compute an inclusion of the derivative of f for intervals
  x. f and df must be functions of one scalar variable, or
  expressions in x. Additional parameters after the semicolon are
  passed to both functions.
  
  The initial interval x must already contain a zero. If x is a
  point, and not an interval, the function tries to get an initial
  interval with the usual Newton method.
  
  Returns {x0,f}: the solution and a flag, if the solution is
  guaranteed.
  
  >expr &= x^2*cos(x);
  >inewton(expr,&diff(expr,x),1,y=0.1)
  ~-0.32483576255267282,-0.32483576255267244~
  >longest solve(expr,0.3,y=0.1)
  0.3248357625526727
  
  See: 
inewton2 (Numerical Algorithms),
mxminewton (Maxima Functions for Euler),
inewton2 (Numerical Algorithms)
function inewton2 (f$:string, Df$:string, x:interval, check:integer=false)
  Guaranteed inclusion of the zero of f, a function of several parameters.
  
  Works like newton2, starting from a interval vector x which already
  contains a solution. If x is no interval, the function tries to
  find such an interval.
  
  f and Df must be a function of a row vector x, or an expression in
  x. f must return a row vector, and Df the Jacobi matrix of f.
  
  Returns {x,valid}.
  
  If check is false, the result is not checked for a guaranteed
  inclusion. In this case the return value of valid can be checked
  to learn, if the inclusion is a guaranteed inclusion. If checked
  is true valid=0 will throw an error exception.
  
function plotintervals (r, style=none)
  Adds plots of two dimensional intervals to a given plot.
  
  r is an nx2 vector of intervals containing the x values in the
  first row, and the y values in the second row. Each value is
  plotted as a bar.
  
  plot2d() can also plot (x,y), where y is a vector of intervals. But
  there the values of y are assumed to be the result of a function
  with errors.
  
  >expr &= x^3-x;
  >plot2d(expr,r=1.1);
  >x=(-1:0.1:1)' + ~-0.02,0.02~;
  >plotintervals(x|expr(x),style="O"):
  
  See: 
mxmibisectfxy (Maxima Functions for Euler)
function arnoldi (A, q ,m)
  Compute an orthonormal basis of the Krylov space
  
  This function computes a matrix Q with orthonormal
  columns which are a basis of the Krylov matrix
  
  b, Ab, A^2.b, ..., A^m.b
  
  A : real nxk matrix
  q : real nx1 vector
  m : integer
  
  Result Q : real nx(m+1) matrix
  
  n=6; A=normal(n,n); q=normal(n,1); ...
  Q = arnoldi(A,q,4)
  Q'.Q // -> id(5)
  Q'.A.Q // -> Hessenberg matrix
  

Documentation Homepage