Sysquake Pro – Table of Contents

Sysquake for LaTeX – Table of Contents

# 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)``c=0.1``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

#### See also

`optimset`,
`fminsearch`,
`fzero`,
`inline`,
operator `@`

### 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``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

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`

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

#### See also

`optimset`,
`fminbnd`,
`lsqnonlin`,
`fsolve`,
`inline`,
operator `@`

### 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

#### See also

`optimset`,
`fminsearch`,
`fzero`,
`inline`,
operator `@`

### 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)``c=10``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

#### See also

`optimset`,
`fminsearch`,
`inline`,
operator `@`,
`roots`

### 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

`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

#### See also

`sum`,
`ode45`,
`inline`,
operator `@`

### 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_');

#### See also

`optimset`,
`lsqnonlin`,
`inline`,
operator `@`

### 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_');

#### See also

`optimset`,
`fminsearch`,
`lsqcurvefit`,
`inline`,
operator `@`

### 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``y2=x'`

`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``x'(0)=0``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``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)]);

#### See also

`odeset`,
`integral`,
`inline`,
operator `@`,
`expm`

### 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

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

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)``u(t)=0``t<1``u(t)=1``t>=1`

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``u=0``filterfun` is *less than or equal to* and not *less than*.
Then the event occurs, and integration continues from
`t=1+eps``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

(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)``d``t>0``x'(t)=Ax(t)+Bu(t)``y(t)=Cx(t-d)``u(t)=1-y(t)``x'(t)=Ax(t)+B(1-Cx(t-d))`

Delayed state is interpolated from past results with `interp1`. Note that values
for `t<0``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).';

#### See also

`ode23`,
`ode45`,
`ode23s`,
`optimset`,
`interp1`

### 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)`

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

#### See also

`fzero`,
`fminbnd`,
`fminsearch`,
`lsqnonlin`,
`lsqcurvefit`

### 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

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