Non-Linear Numerical Functions

fminbnd

Minimum of a function.

Syntax

(x, y) = fminbnd(fun, x0)
(x, y) = fminbnd(fun, [xlow,xhigh])
(x, y) = fminbnd(..., options)
(x, y) = fminbnd(..., options, ...)
(x, y, didConverge) = fminbnd(...)

Description

fminbnd(fun,...) finds numerically a local minimum of function fun. fun is either specified by its name or given as an anonymous or inline function or a function reference. It has at least one input argument x, and it returns one output argument, also a real number. fminbnd finds the value x such that fun(x) is minimized.

Second argument tells where to search; it can be either a starting point or a pair of values which must bracket the minimum.

The optional third argument may contain options. It is either the empty array [] for default options, or the result of optimset.

Remaining input arguments of fminbnd, if any, are given as additional input arguments to function fun. They permit to parameterize the function. For example fminbnd('fun',x0,[],2,5) calls fun as fun(x,2,5) and minimizes its value with respect to x.

The first output argument of fminbnd is the value of x at optimum. The second output argument, if it exists, is the value of fun(x) at optimum. The third output argument, if it exists, is set to true if fminbnd has converged to an optimum, or to false if it has not; in that case, other output arguments are set to the best value obtained. With one or two output arguments, fminbnd throws an error if it does not converge.

Examples

Minimum of a sine near 2, displayed with 15 digits:

fprintf('%.15g\n', fminbnd(@sin, 2));
  4.712389014989218

To find the minimum of c*exp(x)-sin(x) between -1 and 10 with c=0.1, the expression is written as an inline function stored in variable fun:

fun = inline('c*exp(x)-sin(x)', 'x', 'c');

Then fminbnd is used, with the value of c passed as an additional argument:

x = fminbnd(fun,[-1,10],[],0.1)
  x =
    1.2239

With an anonymous function, this becomes

c = 0.1;
fun = @(x) c*exp(x)-sin(x);
x = fminbnd(fun,[-1,10])
  x =
    1.2239

Attempt to find the minimum of an unbounded function:

(x,y,didConverge) = fminbnd(@exp,10)
  x =
    -inf
  y =
    0
  didConverge =
    false

fminsearch

Minimum of a function in R^n.

Syntax

x = fminsearch(fun, x0)
x = fminsearch(..., options)
x = fminsearch(..., options, ...)
(x, y, didConverge) = fminsearch(...)

Description

fminsearch(fun,x0,...) finds numerically a local minimum of function fun. fun is either specified by its name or given as an anonymous or inline function or a function reference. It has at least one input argument x, a real scalar, vector or array, and it returns one output argument, a scalar real number. fminsearch finds the value x such that fun(x) is minimized, starting from point x0.

The optional third input argument may contain options. It is either the empty array [] for default options, or the result of optimset.

Remaining input arguments of fminsearch, if any, are given as additional input arguments to function fun. They permit to parameterize the function. For example fminsearch('fun',x0,[],2,5) calls fun as fun(x,2,5) and minimizes its value with respect to x.

The first output argument of fminsearch is the value of x at optimum. The second output argument, if it exists, is the value of fun(x) at optimum. The third output argument, if it exists, is set to true if fminsearch has converged to an optimum, or to false if it has not; in that case, other output arguments are set to the best value obtained. With one or two output arguments, fminsearch throws an error if it does not converge.

Algorithm

fminsearch implements the Nelder-Mead simplex method. It starts from a polyhedron centered around x0 (the "simplex"). Then at each iteration, either vertex x_i with the maximum value fun(x_i) is moved to decrease it with a reflexion-expansion, a reflexion, or a contraction; or the simplex is shrinked around the vertex with minimum value. Iterations stop when the simplex is smaller than the tolerance, or when the maximum number of iterations or function evaluations is reached (then an error is thrown).

Examples

Minimum of a sine near 2, displayed with 15 digits:

fprintf('%.15g\n', fminsearch(@sin, 2));
  4.712388977408411

Maximum of x*exp(-(x*y)^2)*x*y-0.1*x^2 The function is defined as an anonymous function stored in variable fun:

fun = @(x,y) x.*exp(-(x.*y).^2).*x.*y-0.1*x.^2;

In Sysquake, the contour plot can be displayed with the following commands:

[X,Y] = meshgrid(0:0.02:3, 0:0.02:3);
contour(feval(fun, X, Y), [0,3,0,3], 0.1:0.05:0.5);

The maximum is obtained by minimizing the opposite of the function, rewritten to use as input a single variable in R^2:

mfun = @(X) -(X(1)*exp(-(X(1)*X(2))^2)*X(1)*X(2)-0.1*X(1)^2);
fminsearch(mfun, [1, 2])
  2.1444  0.3297

Here is another way to find this maximum, by calling fun from an intermediate anonymous function:

fminsearch(@(X) -fun(X(1),X(2)), [1, 2])
  2.1444  0.3297

For the same function with a constraint x<1, the objective function can be modified to return infinity for inputs outside the feasible region (note that we can start on the constraint boundary, but starting from the infeasible region would probably fail):

fminsearch(@(X) X(1) < 1 ? -fun(X(1),X(2)) : inf, [1, 2])
  1    0.7071

fsolve

Solve a system of nonlinear equations.

Syntax

x = fsolve(fun, x0)
x = fsolve(..., options)
x = fsolve(..., options, ...)
(x, y, didConverge) = fsolve(...)

Description

fsolve(fun,x0,...) finds numerically a zero of function fun. fun is either specified by its name or given as an anonymous or inline function or a function reference. It has at least one input argument x, a real scalar, vector or array, and it returns one output argument y whose size should match x. fsolve attempts to find the value x such that fun(x) is zero, starting from point x0. Depending on the existence of any solution and on the choice of x0, fsolve may fail to find a zero.

The optional third input argument may contain options. It is either the empty array [] for default options, or the result of optimset.

Remaining input arguments of fsolve, if any, are given as additional input arguments to function fun. They permit to parameterize the function. For example fsolve(@fun,x0,[],2,5) finds the value of x such that the result of fun(x,2,5) is zero.

The first output argument of fsolve is the value of x at zero. The second output argument, if it exists, is the value of fun(x) at zero; it should be a vector or array whose elements are zero, up to the tolerance, unless fsolve cannot find it. The third output argument, if it exists, is set to true if fsolve has converged to a solution, or to false if it has not; in that case, other output arguments are set to the best value obtained. With one or two output arguments, fsolve throws an error if it does not converge.

Algorithm

fsolve minimizes the sum of squares of the vector elements returned by fun using the Nelder-Mead simplex method of fminsearch.

Example

One of the zeros of x1^2+x2^2=10, x2=exp(x1):

[x, y, didConverge] = fsolve(@(x) [x(1)^2+x(2)^2-10; x(2)-exp(x(1))], [0; 0])
  x =
   -3.1620
    0.0423
  y =
   -0.0000
   -0.0000
  didConverge =
    true

fzero

Zero of a function.

Syntax

x = fzero(fun,x0)
x = fzero(fun,[xlow,xhigh])
x = fzero(...,options)
x = fzero(...,options,...)

Description

fzero(fun,...) finds numerically a zero of function fun. fun is either specified by its name or given as an anonymous or inline function or a function reference. It has at least one input argument x, and it returns one output argument, also a real number. fzero finds the value x such that fun(x)==0, up to some tolerance.

Second argument tells where to search; it can be either a starting point or a pair of values xlow and xhigh which must bracket the zero, such that fun(xlow) and fun(xhigh) have opposite sign.

The optional third argument may contain options. It is either the empty array [] for the default options, or the result of optimset.

Additional input arguments of fzero are given as additional input arguments to the function specified by fun. They permit to parameterize the function.

Examples

Zero of a sine near 3, displayed with 15 digits:

fprintf('%.15g\n', fzero(@sin, 3));
  3.141592653589793

To find the solution of exp(x)=c+sqrt(x) between 0 and 100 with c=10, a function f whose zero gives the desired solution is written:

function y = f(x,c)
  y = exp(x) - c - sqrt(x);

Then fsolve is used, with the value of c passed as an additional argument:

x = fzero(@f,[0,100],[],10)
  x =
    2.4479
f(x,10)
  1.9984e-15

An anonymous function can be used to avoid passing 10 as an additional argument, which can be error-prone since a dummy empty option arguments has to be inserted.

x = fzero(@(x) f(x,10), [0,100])
  x =
    2.4479

integral

Numerical integration.

Syntax

y = integral(fun, a, b)
y = integral(fun, a, b, options)

Description

integral(fun,a,b) integrates numerically function fun between a and b. fun is either specified by its name or given as an anonymous or inline function or a function reference. It has a single input argument and a single output argument, both scalar real or complex.

Options can be provided with named arguments. The following options are accepted:

Name	Default	Meaning
`AbsTol`	`1e-6`	maximum absolute error
`RelTol`	`1e-3`	maximum relative error
`Display`	`false`	statistics display

Example

Integration of t*exp(-t) between 0 and 2 with an absolute tolerance of 1e-9:

integral(@(t) t*exp(-t), 0, 2, AbsTol=1e-9)
  0.5940

lsqcurvefit

Least-square curve fitting.

Syntax

param = lsqcurvefit(fun, param0, x, y)
param = lsqcurvefit(..., options)
param = lsqcurvefit(..., options, ...)
(param, r, didConverge) = lsqcurvefit(...)

Description

lsqcurvefit(fun,p0,x,y,...) finds numerically the parameters of function fun such that it provides the best fit for the curve defined by x and y in a least-square sense. fun is either specified by its name or given as an anonymous or inline function or a function reference. It has at least two input arguments: p, the parameter vector, and x, a vector or array of input data; it returns one output argument, a vector or array the same size as x and y. Its header could be

function y = f(param, x)

lsqcurvefit finds the value p which minimizes sum((fun(p,x)-y).^2), starting from parameters p0. All values are real.

The optional fifth input argument may contain options. It is either the empty array [] for default options, or the result of optimset.

Remaining input arguments of lsqcurvefit, if any, are given as additional input arguments to function fun. They permit to parameterize the function. For example lsqcurvefit('fun',p0,x,y,[],2,5) calls fun as fun(p,x,2,5) and find the (local) least-square solution with respect to p.

The first output argument of lsqcurvefit is the value of p at optimum. The second output argument, if it exists, is the value of the cost function at optimum. The third output argument, if it exists, is set to true if lsqcurvefit has converged to an optimum, or to false if it has not; in that case, other output arguments are set to the best value obtained. With one or two output arguments, lsqcurvefit throws an error if it does not converge.

Algorithm

Like lsqnonlin, lsqcurvefit is based on the Nelder-Mead simplex method.

Example

Find the best curve fit of y=a*sin(b*x+c) with respect to parameters a, b and c, where x and y are given (see the example of lsqnonlin for another way to solve the same problem).

% assume nominal parameter values a0=2, b0=3, c0=1
a0 = 2; b0 = 3; c0 = 1;
% reset the seed of rand and randn for reproducible results
rand('s', 0); randn('s', 0);
% create x and y, with noise
x0 = rand(1, 100);
x = x0 + 0.05 * randn(1, 100);
y = a0 * sin(b0 * x0 + c0) + 0.05 * randn(1, 100);
% find least-square curve fit, starting from 1, 1, 1
p0 = [1; 1; 1];
p_ls = lsqcurvefit(@(p, x) p(1) * sin(p(2) * x + p(3)), p0, x, y)
  p_ls =
    2.0060
    2.8504
    1.0836

In Sysquake, the solution can be displayed with

fplot(@(x) a0 * sin(b0 * x + c0), [0,1], 'r');
plot(x, y, 'o');
fplot(@(x) p_ls(1)*sin(p_ls(2)*x+p_ls(3)), [min(x), max(x)]);
legend('Nominal\nSamples\nLS fit', 'r_kok_');

lsqnonlin

Nonlinear least-square solver.

Syntax

x = lsqnonlin(fun, x0)
x = lsqnonlin(..., options)
x = lsqnonlin(..., options, ...)
(x, y, didConverge) = lsqnonlin(...)

Description

lsqnonlin(fun,x0,...) finds numerically the value such that the sum of squares of the output vector produced by fun is a local minimum. fun is either specified by its name or given as an anonymous or inline function or a function reference. It has at least one input argument x, a real scalar, vector or array, and it returns one output argument, a real vector or array. Its header could be

function y = f(x)

lsqnonlin finds the value x such that sum(fun(x(:)).^2) is minimized, starting from point x0.

The optional third input argument may contain options. It is either the empty array [] for default options, or the result of optimset.

Remaining input arguments of lsqnonlin, if any, are given as additional input arguments to function fun. They permit to parameterize the function. For example lsqnonlin('fun',x0,[],2,5) calls fun as fun(x,2,5) and find the (local) least-square solution with respect to x.

The first output argument of lsqnonlin is the value of x at optimum. The second output argument, if it exists, is the value of fun(x) at optimum. The third output argument, if it exists, is set to true if lsqnonlin has converged to an optimum, or to false if it has not; in that case, other output arguments are set to the best value obtained. With one or two output arguments, lsqnonlin throws an error if it does not converge.

Algorithm

Like fminsearch, lsqnonlin is based on the Nelder-Mead simplex method.

Example

Find the least-square solution of a*sin(b*x+c)-y with respect to parameters a, b and c, where x and y are given (see the example of lsqcurvefit for another way to solve the same problem).

% assume nominal parameter values a0=2, b0=3, c0=1
a0 = 2; b0 = 3; c0 = 1;
% reset the seed of rand and randn for reproducible results
rand('s', 0); randn('s', 0);
% create x and y, with noise
x0 = rand(1, 100);
x = x0 + 0.05 * randn(1, 100);
y = a0 * sin(b0 * x0 + c0) + 0.05 * randn(1, 100);
% find least-square solution, starting from 1, 1, 1
p0 = [1; 1; 1];
p_ls = lsqnonlin(@(p) p(1) * sin(p(2) * x + p(3)) - y, p0)
  p_ls =
    2.0060
    2.8504
    1.0836

In Sysquake, the solution can be displayed with

fplot(@(x) a0 * sin(b0 * x + c0), [0,1], 'r');
plot(x, y, 'o');
fplot(@(x) p_ls(1)*sin(p_ls(2)*x+p_ls(3)), [min(x), max(x)]);
legend('Nominal\nSamples\nLS fit', 'r_kok_');

ode23 ode45 ode23s

Ordinary differential equation integration.

Syntax

(t,y) = ode23(fun,[t0,tend],y0)
(t,y) = ode23(fun,[t0,tend],y0,options)
(t,y) = ode23(fun,[t0,tend],y0,options,...)
(t,y,te,ye,ie) = ode23(...)
(t,y) = ode45(fun,[t0,tend],y0)
(t,y) = ode23s(fun,[t0,tend],y0)
...

Description

ode23(fun,[t0,tend],y0) and ode45(fun,[t0,tend],y0) integrate numerically an ordinary differential equation (ODE). Both functions are based on a Runge-Kutta algorithm with adaptive time step; ode23 is low-order and ode45 high-order. In most cases for non-stiff equations, ode45 is the best method.

ode23s(fun,[t0,tend],y0) integrates numerically an ordinary differential equation with a low-order algorithm suitable for stiff systems.

The function to be integrated is either specified by its name or given as an anonymous or inline function or a function reference. It should have at least two input arguments and exactly one output argument:

function yp = f(t,y)

The function calculates the derivative yp of the state vector y at time t.

Integration is performed over the time range specified by the second argument [t0,tend], starting from the initial state y0. It may stop before reaching tend if the integration step cannot be reduced enough to obtain the required tolerance. If the function is continuous, you can try to reduce MinStep in the options argument (see below).

The optional fourth argument may contain options. It is either the empty array [] for the default options, or the result of odeset (the use of a vector of option values is deprecated.)

Events generated by options Events or EventTime can be obtained by three additional output arguments: (t,y,te,ye,ie)=... returns event times in te, the corresponding states in ye and the corresponding event identifiers in ie.

Additional input arguments of ode45 are given as additional input arguments to the function specified by fun. They permit to parameterize the ODE.

ode23s needs the jacobian of the ODE. The jacobian can be passed as a constant square matrix or as a function in the Jacobian option; otherwise numerical approximations are computed.

Examples

Let us integrate the following ordinary differential equation (Van Der Pol equation), parameterized by mu:

x" = mu (1 - x^2) x' - x

Let y1=x and y2=x'; their derivatives are

y1' = y2
y2' = mu (1 - y1^2) y2 - y1

and can be computed by the following function:

function yp = f(t, y, mu)
yp = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

The following ode45 call integrates the Van Der Pol equation from 0 to 10 with the default options, starting from x(0)=2 and x'(0)=0, with mu=1:

(t, y) = ode45(@f, [0,10], [2;0], [], 1);

The same result can be obtained with an anonymous function:

mu=1;
(t, y) = ode45(@(t,y) [y(2); mu*(1-y(1)^2)*y(2)-y(1)],
  [0,10], [2;0]);

The plot command expects traces along the second dimension; consequently, the result of ode45 should be transposed.

plot(t', y');

If mu is large, the ODE is stiff. ode23s is much more efficient.

mu=100;
(t, y) = ode23s(@(t,y) [y(2); mu*(1-y(1)^2)*y(2)-y(1)],
  [0,300], [2;0]);

While ode23s can be called with the same arguments as ode23 and ode45, it is more efficient to provide a function which computes directly the jacobian of the ODE to avoid numerical approximations. The jacobian of the Van Der Pol equations can be computed by an anonymous function and passed to ode23s as a named argument:

mu=100;
(t, y) = ode23s(@(t,y) [y(2); mu*(1-y(1)^2)*y(2)-y(1)],
  [0,300], [2;0],
  Jacobian=@(t,y) [0, 1; -2*mu*y(1)*y(2)-1, mu*(1-y(1)^2)]);

odeset

Options for ordinary differential equation integration.

Syntax

options = odeset
options = odeset(name1=value1, ...)
options = odeset(name1, value1, ...)
options = odeset(options0, name1=value1, ...)
options = odeset(options0, name1, value1, ...)

Description

odeset(name1,value1,...) creates the option argument used by ode23, ode45 and ode23s. Options are specified with name/value pairs, where the name is a string which must match exactly the names in the table below. Case is significant. Alternatively, options can be given with named arguments. Options which are not specified have a default value. The result is a structure whose fields correspond to each option. Without any input argument, odeset creates a structure with all the default options. Note that ode23 etc. also interpret the lack of an option argument, or the empty array [], as a request to use the default values. Options can also be passed directly to them as named arguments.

When its first input argument is a structure, odeset adds or changes fields which correspond to the name/value pairs which follow.

Here is the list of permissible options (empty arrays mean "automatic"):

Name	Default	Meaning
`AbsTol`	`1e-6`	maximum absolute error
`Events`	`[]` (none)	state-based event function
`EventTime`	`[]` (none)	time-based event function
`InitialStep`	`[]` (`10*MinStep`)	initial time step
`Jacobian`	`[]` (undefined)	ODE jacobian
`MaxStep`	`[]` (time range/10)	maximum time step
`MinStep`	`[]` (time range/1e6)	minimum time step
`NormControl`	`false`	error control on state norm
`OnEvent`	`[]` (none)	event function
`OutputFcn`	`[]` (none)	output function
`Past`	`false`	provide past times and states
`PreArg`	`{}`	list of prepended input arguments
`Refine`	`[]` (1, 4 for `ode45`)	refinement factor
`RelTol`	`1e-3`	maximum relative error
`Stats`	`false`	statistics display

Jacobian

The jacobian is used only by ode23s. With an empty matrix, ode23s calculates numerical approximations. If possible, it is more efficient to provide the jacobian either as a constant square matrix if it is constant, or as a function called by ode23s at each integration step. The function is defined as

function J = jac(t, y)

The jacobian of vector function f(y) is the square matrix whose columns are the partial derivatives of f with respect to y(1), y(2) etc.

Time steps and output

Several options control how the time step is tuned during the numeric integration. Error is calculated separately on each element of y if NormControl is false, or on norm(y) if it is true; time steps are chosen so that it remains under AbsTol or RelTol times the state, whichever is larger. If this cannot be achieved, for instance if the system is stiff and requires an integration step smaller than MinStep, integration is aborted.

'Refine' specifies how many points are added to the result for each integration step. When it is larger than 1, additional points are interpolated, which is much faster than reducing MaxStep.

The output function OutputFcn, if defined, is called after each step. It is a function name in a string, a function reference, or an anonymous or inline function, which can be defined as

function stop = outfun(tn, yn)

where tn is the time of the new samples, yn their values, and stop a logical value which is false to continue integrating or true to stop. The number of new samples is given by the value of Refine; when multiple values are provided, tn is a row vector and yn is a matrix whose columns are the corresponding states. The output function can be used for incremental plots, for animations, or for managing large amounts of output data without storing them in variables.

Events

Events are additional time steps at controlled time, to change instantaneously the states, and to base the termination condition on the states. Time instants where events occur are either given explicitly by EventTime, or implicitly by Events. There can be multiple streams of events, which are checked independently and are identified by a positive integer for Events, or a negative integer for EventTime. For instance, for a ball which bounces between several walls, the intersection between each wall and the ball trajectory would be a different event stream.

For events which occur at regular times, EventTime is an n-by-two matrix: for each row, the first column gives the time step ts, and the second column gives the offset to. Non-repeating events are specified with an infinite time step ts. Events occur at time t=to+k*ts, where k is an integer.

When event time is varying, EventTime is a function which can be defined as

function eventTime = eventtimefun(t, y, ...)

where t is the current time, y the current state, and the ellipsis stand for additional arguments passed to ode*. The function returns a (column) vector whose elements are the times where the next event occurs. In both cases, each row corresponds to a different event stream.

For events which are based on the state, the value of a function which depends on the time and the states is checked; the event occurs when its sign changes. Events is a function which can be defined as

function (value, isterminal, direction) ...
          = eventsfun(t, y, ...)

Input arguments are the same as for EventTime. Output arguments are (column) vectors where each element i corresponds to an event stream. An event occurs when value(i) crosses zero, in either direction if direction(i)==0, from negative to nonnegative if direction(i)>0, or from positive to nonpositive if direction(i)<0. The event terminates integration if isterminal(i) is true. The Events function is evaluated for each state obtained by integration; intermediate time steps obtained by interpolation when Refine is larger than 1 are not considered. When an event occurs, the integration time step is reset to the initial value, and new events are disabled during the next integration step to avoid shattering. MaxStep should be used if events are missed when the result of Events is not monotonous between events.

When an event occurs, function OnEvent is called if it exists. It can be defined as

function yn = onevent(t, y, i, ...)

where i identifies the event stream (positive for events produced by Events or negative for events produced by EventTime); and the output yn is the new value of the state, immediately after the event.

The primary goal of ode* functions is to integrate states. However, there are systems where some states are constant between events, and are changed only when an event occurs. For instance, in a relay with hysteresis, the output is constant except when the input overshoots some value. In the general case, ni states are integrated and n-ni states are kept constant between events. The total number of states n is given by the length of the initial state vector y0, and the number of integrated states ni is given by the size of the output of the integrated function. Function OnEvent can produce a vector of size n to replace all the states, of size n-ni to replace the non-integrated states, or empty to replace no state (this can be used to display results or to store them in a file, for instance).

Event times are computed after an integration step has been accepted. If an event occurs before the end of the integration step, the step is shortened; event information is stored in the output arguments of ode* te, ie and ye; and the OnEvent function is called. The output arguments t and y of ode* contain two rows with the same time and the state right before the event and right after it. The time step used for integration is not modified by events.

Additional arguments

Past is a logical value which, if true, specifies that the time and state values computed until now (what will eventually be the result of ode23, ode45 or ode23s) are passed as additional input arguments to functions called during intergration. This is especially useful for delay differential equations (DDE), where the state at some time point in the past can be interpolated from the integration results accumulated until now with interp1. Assuming no additional parameters or PreArg (see below), functions must be defined as

function yp = f(t,y,tpast,ypast)
function stop = outfun(tn,yn,tpast,ypast)
function eventTime = eventtimefun(t,y,tpast,ypast)
function (value, isterminal, direction) ...
          = eventsfun(t,y,tpast,ypast)
function yn = onevent(t,y,tpast,ypast,i)
function J = jac(t,y,tpast,ypast)

PreArg is a list of additional input arguments for all functions called during integration; they are placed before normal arguments. For example, if its value is {1,'abc'}, the integrated function is called as fun(1,'abc',t,y), the output function as outfun(1,'abc',tn,yn), and so on.

Examples

Default options

odeset
  AbsTol: 1e-6
  Events: []
  EventTime: []
  InitialStep: []
  Jacobian: []
  MaxStep: []
  MinStep: []
  NormControl: false
  OnEvent: []
  OutputFcn: []
  PreArg: {}
  Refine: []
  RelTol: 1e-3
  Stats: false

Options passed as named arguments

Unless options must be stored as a whole in a variable, it is often more convenient to pass them directly to the integration function as named arguments. The following calls are equivalent.

(t, y) = ode45(fun, tspan, y0, odeset('RelTol', 1e-4));
(t, y) = ode45(fun, tspan, y0, odeset(RelTol=1e-4));
(t, y) = ode45(fun, tspan, y0, RelTol=1e-4);

Option `Refine`

ode45 is typically able to use large time steps to achieve the requested tolerance. When plotting the output, however, interpolating it with straight lines produces visual artifacts. This is why ode45 inserts 3 interpolated points for each calculated point, based on the fifth-order approximation calculated for the integration (Refine is 4 by default). In the following code, curves with and without interpolation are compared. Note that the numbers of evaluations of the function being integrated are the same.

mu = 1;
fun = @(t,y) [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
(t, y) = ode45(fun, [0,5], [2;0],
               Refine=1, Stats=true);
  Number of function evaluations: 289
  Successful steps: 42
  Failed steps (error too large): 6
size(y)
  43  2
(ti, yi) = ode45(fun, [0,5], [2;0],
                 Stats=true);
  Number of function evaluations: 289
  Successful steps: 42
  Failed steps (error too large): 6
size(yi)
  169  2
plot(ti', yi', 'g');
plot(t', y');

State-based events

For simulating a ball bouncing on the ground, an event is generated every time the ball hits the ground, and its speed is changed instantaneously. Let y(1) be the height of the ball above the ground, and y(2) its speed (SI units are used). The state-space model is

  y' = [y(2); -9.81];

An event occurs when the ball hits the ground:

value = y(1);
isterminal = false;
direction = -1;

When the event occurs, a new state is computed:

yn = [0; -damping*y(2)];

To integrate this, the following functions are defined:

function yp = ballfun(t, y, damping)
  yp = [y(2); -9.81];
function (v, te, d) = ballevents(t, y, damping)
  v = y(1);    // event when the height becomes negative
  te = false;  // do not terminate
  d = -1;      // only for negative speeds
function yn = ballonevent(t, y, i, damping)
  yn = [0; -damping*y(2)];

Ball state is integrated during 5 s with

opt = odeset(Events=@ballevents,
             OnEvent=@ballonevent);
(t, y) = ode45(@ballfun, [0, 5], [2; 0], opt, 1);
plot(t', y');

Time events with discontinuous function

If the function being integrated has discontinuities at known time instants, option EventTime can be used to insure an accurate switching time. Consider a first-order filter with input u(t), where u(t)=0 for t<1 and u(t)=1 for t>=1. The following function is defined for the state derivative:

function yp = filterfun(t, y)
  yp = -y + (t <= 1 ? 0 : 1);

A single time event is generated at t=1:

opt = odeset(EventTime=[inf, 1]);
(t, y) = ode45(@filterfun, [0, 5], 0, opt);
plot(t', y');

Function filterfun is integrated in the normal way until t=1 inclusive, with u=0. This is why the conditional expression in filterfun is less than or equal to and not less than. Then the event occurs, and integration continues from t=1+eps with u=1.

Early termination

The normal termination criterion is the final time specified in the tspan argument of ODE solver functions. A state-based criterion can be specified with either a state-based event Events or an output function OutputFcn. An output function can be simpler to specify, and faster because it does not attempt to reach precisely a state-based transition.

The next example integrates the free fall of an object until its height becomes negative. The state contains the height in x(1) and the height derivative (the speed) in x(2). The objects starts at rest at an height of 10 m. The final time specified in tspan (100 s) is assumed to be large enough. The integration terminates as soon as any new height is negative (multiple samples are fed to OutputFcn at each integration step if Refine>1).

fder = @(t,x) [x(2); -9.81];
x0 = [10; 0];
(t,x) = ode45(fder, [0,100], x0, OutputFcn=@(tn,xn) any(xn(:,1)<0));

Non-integrated state

For the last example, we will consider a system made of an integrator and a relay with hysteresis in a loop. Let y(1) be the output of the integrator and y(2) the output of the relay. Only y(1) is integrated:

yi' = y(2);

An event occurs when the integrator is larger or smaller than the hysteresis:

value = y(1) - y(2);
isTerminal = false;
direction = sign(y(2));

When the event occurs, a new value is computed for the 2nd state:

yn = -y(2);

To integrate this, the following functions are defined:

function yp = relayfun(t, y)
  yp = y(2);
function (v, te, d) = relayevents(t, y)
  v = y(1) - y(2);
  te = false;
  d = sign(y(2));
function yn = relayonevent(t, y, i)
  yn = -y(2);

The initial state is [0;1]; 0 for the integrator, and 1 for the output of the relay. State is integrated during 5 s with

(t, y) = ode45(@relayfun, [0, 5], [0; 1],
               Events=@relayevents, OnEvent=@relayonevent);
plot(t', y');

Delay differential equation

A system whose Laplace transform is Y(s)/U(s)=exp(-ds)/(s^2+s) (first order + integrator + delay d) is simulated with unit negative feedback. The reference signal is 1 for t>0. First, the open-loop system is converted from transfer function to state-space, such that x'(t)=Ax(t)+Bu(t) and y(t)=Cx(t-d). The closed-loop state-space model is obtained by setting u(t)=1-y(t), which gives x'(t)=Ax(t)+B(1-Cx(t-d)).

Delayed state is interpolated from past results with interp1. Note that values for t<0 (extrapolated) are set to 0, and that values more recent than the last result are interpolated with the state passed to f for current t.

(A,B,C) = tf2ss(1,[1,1,0]);
d = 0.1;
x0 = zeros(length(A),1);
tmax = 10;
f = @(t,x,tpast,xpast) ...
     A*x+B*(1-C*interp1([tpast;t],[xpast;x.'],t-d,'1',0).');
(t,x) = ode45(f, [0,tmax], x0, Past=true);

Output y can be computed from the state:

y = C * interp1(t,x,t-d,'1',0).';

optimset

Options for minimization and zero finding.

Syntax

options = optimset
options = optimset(name1=value1, ...)
options = optimset(name1, value1, ...)
options = optimset(options0, name1=value1, ...)
options = optimset(options0, name1, value1, ...)

Description

optimset(name1,value1,...) creates the option argument used by fminbnd, fminsearch, fzero, fsolve, and other optimization functions. Options are specified with name/value pairs, where the name is a string which must match exactly the names in the table below. Case is significant. Alternatively, options can be given with named arguments. Options which are not specified have a default value. The result is a structure whose fields correspond to each option. Without any input argument, optimset creates a structure with all the default options. Note that fminbnd, fminsearch, and fzero also interpret the lack of an option argument, or the empty array [], as a request to use the default values. Options can also be passed directly to fminbnd and other similar functions as named arguments.

When its first input argument is a structure, optimset adds or changes fields which correspond to the name/value pairs which follow.

Here is the list of permissible options (empty arrays mean "automatic"):

Name	Default	Meaning
`Display`	`false`	detailed display
`MaxFunEvals`	`1000`	maximum number of evaluations
`MaxIter`	`500`	maximum number of iterations
`TolX`	`[]`	maximum relative error

The default value of TolX is eps for fzero and sqrt(eps) for fminbnd and fminsearch.

Examples

Default options:

optimset
  Display: false
  MaxFunEvals: 1000
  MaxIter: 500
  TolX: []

Display of the steps performed to find the zero of cos(x) between 1 and 2:

fzero(@cos, [1,2], optimset('Display',true))
  Checking lower bound
  Checking upper bound
  Inverse quadratic interpolation 2,1.5649,1
  Inverse quadratic interpolation 1.5649,1.571,2
  Inverse quadratic interpolation 1.571,1.5708,1.5649
  Inverse quadratic interpolation 1.5708,1.5708,1.571
  Inverse quadratic interpolation 1.5708,1.5708,1.571
ans =
  1.5708

quad

Numerical integration.

Syntax

y = quad(fun, a, b)
y = quad(fun, a, b, tol)
y = quad(fun, a, b, tol, trace)
y = quad(fun, a, b, tol, trace, ...)

Description

quad(fun,a,b) integrates numerically real function fun between a and b. fun is either specified by its name or given as an anonymous or inline function or a function reference.

The optional fourth argument is the requested relative tolerance of the result. It is either a positive real scalar number or the empty matrix (or missing argument) for the default value, which is sqrt(eps). The optional fifth argument, if true or nonzero, makes quad displays information at each step.

Additional input arguments of quad are given as additional input arguments to function fun. They permit to parameterize the function.

Example

Integration of t*exp(-t) between 0 and 2:

quad(@(t) t*exp(-t), 0, 2)
  0.5940

Remark

Function quad is obsolete and should be replaced with integral, which supports named options and complex numbers.

Non-Linear Numerical Functions

fminbnd

Syntax

Description

Examples

See also

fminsearch

Syntax

Description

Algorithm

Examples

See also

fsolve

Syntax

Description

Algorithm

Example

See also

fzero

Syntax

Description

Examples

See also

integral

Syntax

Description

Example

See also

lsqcurvefit

Syntax

Description

Algorithm

Example

See also

lsqnonlin

Syntax

Description

Algorithm

Example

See also

ode23 ode45 ode23s

Syntax

Description

Examples

See also

odeset

Syntax

Description

Jacobian

Time steps and output

Events

Additional arguments

Examples

Default options

Options passed as named arguments

Option Refine

State-based events

Time events with discontinuous function

Early termination

Non-integrated state

Delay differential equation

See also

optimset

Syntax

Description

Examples

See also

quad

Syntax

Description

Example

Remark

See also

Option `Refine`