Advanced Solver Syntax
The complete syntax of the ODE solver is as follows and is summarized in
[t, y) = ode23(‘ydot’, tspan, yO, options, pl . p2, … )
where the options argument is created with the new odeset function, and pl , p2, … are optional parameters that can be passed to the function file y dot every time it is called. If these optional parameters are used, but no opt ions are set, use opt ions = [ 1 as a placeholder. The function file ydot can now take additional input arguments. It has the form ydot (t, y, flag, pl , p2, … ), where flag is an argument that notifies the function ydot that the solver is expecting a specific kind of information. We will see an example of the f lag argument shortly. The odeset Function
The odeset function creates an options structure to be supplied to the solver. Its syntax is options = odeset(‘namel’, ‘valuel’ ‘name2’, ‘value2
where name is the name of a property and val ue is the value to be assigned to the property.
A simple example will clarify things. The Re fine property is used to increase the number of output points from the solver by an integer factor II. The default value of n is 1 for all the solvers except ode4 5, whose default value is 4 because of the solver’s large step sizes. Suppose we want to solve the equation
dy . – =smt dt for 0 ::: t :::47l’ with y(O) = O. Define the following function file: function ydot = sinefn(t,y)
ydot = sin(t); Then use the odeset function to set the value of Refine to n = 8, and call the ode45 solver, as shown here:
options = odeset(‘Refine’ ,8); [t, y] = ode45(‘sinefn’, [0, 4*pi], 0, options)’; Many properties can be set with the odeset function. To see a list of these. type odeset.
Another property is the Events property, which has two possible value: on and 0f f. It can be used to locate transitions to, from, or through zeros of a user-defined function. This property can be used to detect when the ODE solution reaches a certain value. To use this property, you must write the ODE file ydot;
to have three outputs, as follows:
[value, isterminal, direction] = ydot(t,y,flag)
The vector value should be programmed to contain the vector describing the event when the ftag equals ‘ events’. Thevectorvalueshouldbeprogrammed to contain the derivatives when the flag is not set to ‘ even t s ‘ or is empty.
val ue is evaluated at the beginning and end of each step, and if any of its elements make transitions to, from, or through zero, the solver determines the time when the transition occurred. The direction of the transition is specified in the vector direct ion. A 1indicates positive direction: a -1indicates negative
_ direction, and a a indicates “don’t care.” The isterminal vector is a logical vector of Is and as that specifies whether a zero crossing of the corresponding element in val ue is a terminal
event. A 1corresponds to a terminal event and halts the solver; a a corresponds to a nonterminal event. For example, suppose we want to simulate a dropped ball bouncing up from the floor. The equation of motion of the ball in free flight is my = -11Ig
or y = -g, where y is the ball’s height above the floor. Put this into state-variable form by defining y, = y and )’2 = y. The equations are
Define the following equation file to use with the solver. We use g = 9.81 rnIs2 and an initial height of 10m. The va 1ue vector here consists of the state vector. The “event” here is a bounce, which occurs when the height is O. function [value, isterminal, direction] = ballode(t, y, flag) if (nargin<3) lisempty(flag)
value = [y(2); -9.81]; elseif flag == ‘events’ % Retu~ns height and velo~ity
value = y;
%Instructs to stop when first \
\variable (height) is zero
isterminal = [1; 0];
% Instructs to detect an event only when the first
% variable (height) is decreasing,
% regardless of velocity
direction = [-1; 0];
error([‘unknown flag”‘flag'” .’]);
We assume that the ball loses 10 percent of its speed each time it bounces. The script file to find the ball’s motion up to the second bounce follows. It calls the solver before and after the bounce.
options = odeset(‘Events’, ‘on’)
[tl, y1] = ode45(‘ballode’, [0,10), [10,0], options);
[t2, y2] = ode45(‘ballode’, [t1(length(t1)), 10], …
t = [tl; t2];y = [y1; y2];
plot (t,Y (:,1)),xlabel (,Time (s)’), …
ylabel(‘Ball Height (m)’),axis([O 5010))
The resulting plot of height versus time is shown in Figure 8.8-1.
Stiff Differential Equations A stiff differential equation is one whose response changes rapidly over a time scale that is short compared to the time scale over which we are interested in the solution. Stiff equations present a challenge to-solve numerically. A small
step size is needed to solve for the rapid changes, but many steps are needed to obtain the solution over the longer time interval, and thus a large error might ,
accumulate. For example, the following equation might be considered “stiff”:
y + 100ly + 1000y = 0 (8.8-1)
The characteristic roots are s = -1 and s = -1000. If the initial conditions are y(O) = I and .Y(O) = 0, the closed-form solution is
y(t) = – (1000e-t _ e-1OOOt)
That part of the response due to the term e-1OOOt disappears after approximately t = 4/1000 = 0.004, but the term due to e:’ does not disappear until after t = 4. Thus it would be difficult for a plot to show the solution accurately, and a numerical solver would need a small step size to compute the rapid changes due. to the e-1OOOt term, much smaller than that required to compute the longer-term
response due to the e:’ term. The result can be a large accumulated error because of the small step size combined with the large number of steps required to obtain the full solution. All of the solvers in MATLAB do rather well with moderately stiff equations,
but if you have trouble solving an equation with one of these solvers, try one of the four solvers specifically designed to handle stiff equations. These are ode15s, a variable-order method; ode23s, a low-order method; ode23tb another low-order method; and ode23 t, a trapezoidal method. Their syntax is identical to the ode23 and ode4 5 solvers.
We have not covered all of the new solver capabilities provided in MATLAB. For more information and additional examples, consult online help.