Minimum of a function.
(x, y) = fminbnd(fun, x0) (x, y) = fminbnd(fun, [xlow,xhigh]) (x, y) = fminbnd(..., options) (x, y) = fminbnd(..., options, ...) (x, y, didConverge) = fminbnd(...)
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.
Minimum of a sine near 2, displayed with 15 digits:
fprintf('%.15g\n', fminbnd(@sin, 2)); 4.712389014989218
To find the minimum of
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
optimset, fminsearch, fzero, inline, operator @
Minimum of a function in R^n.
x = fminsearch(fun, x0) x = fminsearch(..., options) x = fminsearch(..., options, ...) (x, y, didConverge) = fminsearch(...)
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.
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).
Minimum of a sine near 2, displayed with 15 digits:
fprintf('%.15g\n', fminsearch(@sin, 2)); 4.712388977408411
Maximum of
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
For the same function with a constraint
mfunc = @(X) ... X(1) < 1 ... ? -(X(1)*exp(-(X(1)*X(2))^2)*X(1)*X(2) - 0.1*X(1)^2) ... : inf; fminsearch(mfunc, [1, 2]) 1 0.7071
optimset, fminbnd, fzero, inline, operator @
Zero of a function.
x = fzero(fun,x0) x = fzero(fun,[xlow,xhigh]) x = fzero(...,options) x = fzero(...,options,...)
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.
Zero of a sine near 3, displayed with 15 digits:
fprintf('%.15g\n', fzero(@sin, 3)); 3.141592653589793
To find the solution of
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
optimset, fminsearch, inline, operator @, roots
Ordinary differential equation integration.
(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) = ode45(fun,[t0,tend],y0) (t,y) = ode45(fun,[t0,tend],y0,options) (t,y) = ode45(fun,[t0,tend],y0,options,...)
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. 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.)
Additional input arguments of ode45 are given as additional input arguments to the function specified by fun. They permit to parameterize the ODE.
Let us integrate the following ordinary differential equation (Van Der Pol equation),
parameterized by
Let
y2' = mu (1 - y1^2) y2 - y1
and may 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
(t,y)=ode45(@f,[0,10],[2;0],[],1);
The plot command expects traces along the second dimension; consequently, the result of ode45 should be transposed.
plot(t', y');
odeset, quad, inline, operator @, expm
Options for ordinary differential equation integration.
options = odeset options = odeset(name1, value1, ...) options = odeset(options0, name1, value1, ...)
odeset(name1,value1,...) creates the option argument used by ode23 and ode45. 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. 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 and ode45 also interpret the lack of an option argument, or the empty array [], as a request to use the default values.
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 |
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 |
PreArg | {} | list of prepended input arguments |
Refine | [] (1, 4 for ode45) | refinement factor |
RelTol | 1e-3 | maximum relative error |
Stats | false | statistics display |
Several options control how the time step is tuned during the numerical 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 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.
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 with fun(1,'abc',t,y), the output function as outfun(1,'abc',tn,yn), and so on.
odeset AbsTol: 1e-6 Events: [] EventTime: [] InitialStep: [] MaxStep: [] MinStep: [] NormControl: false OnEvent: [] OutputFcn: [] PreArg: {} Refine: [] RelTol: 1e-3 Stats: false
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], ... odeset('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], ... odeset('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');
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');
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
function yp = filterfun(t, y) yp = -y + (t <= 1 ? 0 : 1);
A single time event is generated at
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
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], ... odeset('Events', @relayevents, 'OnEvent', @relayonevent)); plot(t', y');
Options for minimization and zero finding.
options = optimset options = optimset(name1, value1, ...) options = optimset(options0, name1, value1, ...)
optimset(name1,value1,...) creates the option argument used by fminbnd, fminsearch, and fzero. 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. 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.
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.
Default options:
optimset Display: false MaxFunEvals: 1000 MaxIter: 500 TolX: []
Display of the steps performed to find the zero of
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
fzero, fminbnd, fminsearch, odeset
Numerical integration.
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, ...)
quad(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.
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.
quad(@(t) t*exp(-t), 0, 2) 0.5940
sum, ode45, inline, operator @