ode
Description
An ode
object defines a system of ordinary differential equations
or differential algebraic equations to solve.
You can solve initial value problems of the form $$y\text{'}=f\left(t,y\right)$$ or problems that involve a mass matrix, $$M\left(t,y\right)y\text{'}=f\left(t,y\right)$$.
Define aspects of the problem using properties of the ode
object, such as
ODEFcn
, InitialTime
, and
InitialValue
. You can select a specific solver to use, or let
MATLAB^{®} choose an appropriate solver based on properties of the equations. After you
create an ode
object, you can solve the equations using the
solve
or solutionFcn
object functions.
Creation
Description
creates an F
= odeode
object
with default properties.
specifies one or more property values using namevalue arguments. For example, you can
specify the equations to be solved, the initial time for integration, and the value of the
solution at the initial time using the F
= ode(Name=Value
)ODEFcn
,
InitialTime
, and InitialValue
properties.
Properties
Problem Definition
ODEFcn
— Equations to solve
function handle
Equations to solve, specified as a function handle that defines the system of differential equations to be integrated. The function handle can be an anonymous function or a handle to a named function file.
The function dydt = ODEFcn(t,y)
, for a scalar
t
and a column vector y
, must return a column
vector dydt
that corresponds to $$f\left(t,y\right)$$. ODEFcn
must accept at least two input
arguments, t
and y
, even if one of the arguments
is not used in the function. To supply parameter values, define
ODEFcn
as a function that accepts three inputs, dydt =
ODEFcn(t,y,p)
, and then use the Parameters
property
to store parameter values (such as a struct
or
cell
array).
For example, to solve a single equation such as $$y\text{'}=5y3$$, you can use the anonymous function @(t,y)
5*y3
.
For a system of equations, the output of ODEFcn
is a vector.
Each element in the vector is the computed value of the right side of one equation.
For example, consider this system of two equations.
$$\begin{array}{l}y{\text{'}}_{1}={y}_{1}+2{y}_{2}\\ y{\text{'}}_{2}=3{y}_{1}+2{y}_{2}\end{array}$$
An anonymous function that defines these equations is @(t,y)
[y(1)+2*y(2); 3*y(1)+2*y(2)]
.
Example:
F.ODEFcn = @myFcn
specifies a handle to a named function file
myFcn.m
containing the equations.
Example: F.ODEFcn = @(t,y) [y(2); y(1)]
specifies an anonymous
function that defines a system of two equations.
Example: F.ODEFcn = @(t,y,p) 5*y*p(1)3*p(2)
specifies an
anonymous function for a single equation that uses two parameters.
Data Types: function_handle
InitialTime
— Initial time for integration
0
(default)  real scalar
Initial time for integration, specified as a real scalar. The value of
InitialTime
is the beginning of the integration interval where
the initial conditions specified in InitialValue
are applied by
the solver before beginning integration steps.
Example:
F.InitialTime = 10;
Data Types: single
 double
InitialValue
— Value of solution at InitialTime
scalar  vector
Value of solution at InitialTime
, specified as a scalar or
vector. InitialValue
must be a vector with the same length as the
output of ODEFcn
so that an initial value is specified for each
equation defined in ODEFcn
. The value of
InitialTime
is the beginning of the integration interval where
the initial conditions specified in InitialValue
are applied by
the solver before beginning integration steps.
Example:
F.InitialValue = 1;
Example:
F.InitialValue = [1 2 3];
Data Types: single
 double
Parameters
— Equation parameters
array
Equation parameters, specified as an array of any size or data type. The values
you store in Parameters
can be supplied to any of the functions
used for the ODEFcn
, Jacobian
,
MassMatrix
, or EventDefinition
properties
by specifying an extra input argument in the function.
For instance, you can supply parameter values stored in the
Parameters
property to ODEFcn
by
specifying ODEFcn
as a function handle that accepts three inputs,
dydt = ODEFcn(t,y,p)
. For example, F.ODEFcn = @(t,y,p)
5*y*p(1)3*p(2)
specifies an anonymous function for a single equation that
uses two parameters. If you then specify F.Parameters = [2 3]
, then
ODEFcn
uses the parameter values p(1) = 2
and p(2) = 3
whenever the solver calls the function.
Example:
F.Parameters = [0.1 0.5]
specifies two parameter
values.
Jacobian
— Jacobian matrix
odeJacobian
object  matrix  function handle
Jacobian matrix, specified as an odeJacobian
object, matrix, or handle to a function that evaluates the Jacobian. The Jacobian is a
matrix of partial derivatives of the functions that define the system of differential
equations.
$$J=\frac{\partial f}{\partial y}=\left[\begin{array}{ccc}\frac{\partial {f}_{1}}{\partial {y}_{1}}& \frac{\partial {f}_{1}}{\partial {y}_{2}}& \cdots \\ \frac{\partial {f}_{2}}{\partial {y}_{1}}& \frac{\partial {f}_{2}}{\partial {y}_{2}}& \cdots \\ \vdots & \vdots & \end{array}\right]$$
For stiff ODE problems, providing information about the Jacobian matrix is critical for reliability and efficiency of the solver. If you do not provide the Jacobian, then the ODE solver approximates it numerically using finite differences.
For large systems of equations where it is not feasible to provide the entire analytic Jacobian, you can specify the sparsity pattern of the Jacobian matrix instead. The solver uses the sparsity pattern to calculate a sparse Jacobian.
You can specify the value of the Jacobian
property as:
An
odeJacobian
object, which can represent either the Jacobian matrix or its sparsity pattern.A constant matrix with calculated values for $$\frac{\partial f}{\partial y}$$.
A handle to a function that computes the matrix elements and that accepts two input arguments,
dfdy = Fjac(t,y)
. To give the function access to parameter values in theParameters
property, specify a third input argument in the function definition,dfdy = Fjac(t,y,p)
.
If you specify a matrix or function handle, then MATLAB converts it to an odeJacobian
object.
Example: F.Jacobian = @Fjac
specifies the function
Fjac
that evaluates the Jacobian matrix.
Example: F.Jacobian = [0 1; 2 1]
specifies a constant Jacobian
matrix.
Example: F.Jacobian = odeJacobian(SparsityPattern=S)
specifies
the Jacobian sparsity pattern using sparse matrix S
.
EventDefinition
— Events to detect
odeEvent
object
Events to detect, specified as an odeEvent
object. Create an odeEvent
object to define expressions that trigger
an event when they cross zero. You can specify the direction of the zero crossing and
what to do when an event triggers, including the use of a callback function.
MassMatrix
— Mass matrix
odeMassMatrix
object  matrix  function handle
Mass matrix, specified as an odeMassMatrix
object, matrix, or handle to a function that evaluates the
mass matrix.
ode
objects can represent problems of the form $$M\left(t,y\right)\text{\hspace{0.17em}}y\text{'}=f\left(t,y\right)$$, where $$M\left(t,y\right)$$ is a mass matrix that can be full or sparse. The mass matrix encodes
linear combinations of derivatives on the left side of the equation.
When the mass matrix is nonsingular, the equation simplifies to $$y\text{'}={M}^{1}\text{\hspace{0.17em}}f\left(t,y\right)$$ and the ODE has a solution for any initial value. However, it is often more convenient and natural to express the model in terms of the mass matrix directly using $$M\left(t,y\right)\text{\hspace{0.17em}}y\text{'}=f\left(t,y\right)$$, and avoiding the computation of the matrix inverse reduces the storage and execution time needed to solve the problem.
When the mass matrix is singular, then the problem is a system of differential algebraic equations (DAEs). A DAE has a solution only when the initial values are consistent; that is, you must specify the initial slope $$y{\text{'}}_{0}$$ using the
InitialSlope
property such that the initial conditions are all consistent, $$M\left({t}_{0},{y}_{0}\right)y{\text{'}}_{0}=f\left({t}_{0},{y}_{0}\right)$$. If the specified initial conditions are not consistent with theInitialTime
,InitialValue
, andMassMatrix
properties, then the solver treats them as guesses and attempts to compute consistent values for the initial slopes that are close to the guesses before continuing to solve the problem. For more information, see Solve Differential Algebraic Equations (DAEs).
You can specify the value of the MassMatrix
property as:
An
odeMassMatrix
object, which represents the mass matrix and its associated properties. You can specify whether the mass matrix is singular or has state dependence.A constant matrix with calculated values for $$M\left(t,y\right)$$.
A handle to a function that computes the matrix elements and that accepts two input arguments,
M = Mass(t,y)
. To give the function access to parameter values stored in theParameters
property, specify a third input argument in the function definition,M = Mass(t,y,p)
.
If you specify a matrix or function handle, then MATLAB converts it to an odeMassMatrix
object.
Example: F.MassMatrix = @Mass
specifies the function
Mass
that evaluates the mass matrix.
Example: F.MassMatrix = [1 0; 2 1]
specifies a constant mass
matrix.
Example: F.MassMatrix =
odeMassMatrix(MassMatrix=@Mass,StateDependence="strong",SparsityPattern=S)
specifies a statedependent mass matrix and the sparsity pattern when the mass matrix
multiplies a vector.
NonNegativeVariables
— Nonnegative solution components
scalar index  vector of indices
Nonnegative solution components, specified as a scalar index or vector of
indices. Use the NonNegativeVariables
property to specify which
solutions the solver must keep nonnegative. If dydt = ODEFcn(t,y)
,
then the indices correspond to elements in the vector y
. For
example, if you specify a value of 3
, then the solver keeps the
solution component y(3)
nonnegative.
Note
NonNegativeVariables
is not available for the
ode23s
solver. Additionally, for ode15s
,
ode23t
, and ode23tb
the property is not
available for problems that have a mass matrix.
Example: F.NonNegativeVariables = [1 3]
specifies that the first
and third solution components must be kept nonnegative.
Data Types: single
 double
 int8
 int16
 int32
 int64
 uint8
 uint16
 uint32
 uint64
InitialSlope
— Initial value of dy/dt
vector
Initial value of dy/dt
, specified as a vector. Use this
property with the ode15s
and ode23t
solvers
when solving DAEs. The specified vector is the initial slope $$y{\text{'}}_{0}$$ such that $$M\left({t}_{0},{y}_{0}\right)y{\text{'}}_{0}=f\left({t}_{0},{y}_{0}\right)$$. If the specified initial conditions are not consistent with the
InitialTime
, InitialValue
, and
MassMatrix
properties, then the solver treats them as guesses
and attempts to compute consistent values for the initial slopes that are close to the
guesses before continuing to solve the problem. For more information, see Solve Differential Algebraic Equations (DAEs).
Data Types: single
 double
Sensitivity
— Sensitivity analysis information
odeSensitivity
object
Sensitivity analysis information, specified as an odeSensitivity
object. The ODE solver performs sensitivity analysis only when the
Sensitivity
property is nonempty. The analysis is with respect
to the Parameters
property, which must be nonempty and numeric.
See odeSensitivity
for more information.
Example: F.Sensitivity = odeSensitivity
performs sensitivity
analysis on all parameters.
Example: F.Sensitivity = odeSensitivity(ParameterIndices=1:2)
performs sensitivity analysis on the first two parameters.
Solver Control
Solver
— Solver method
"auto"
(default)  "stiff"
 "nonstiff"
 name of a solver
Solver method, specified as one of the values in this table.
Value  Description 

 Automatically selects a solver based on available problem
information and tolerances. The 
 Similar to 
 Similar to 
 RungeKutta (2,3) pair. See 
 RungeKutta (4,5) pair. See 
 RungeKutta 8(7) pair with a 7thorder continuous extension. See

 RungeKutta 9(8) pair with an 8thorder continuous extension. See

 Variablestep, variableorder (VSVO) solver of orders 1 to 13. See

 Variablestep, variableorder (VSVO) solver based on the numerical
differentiation formulas (NDFs) of orders 1 to 5. See 
 Modified Rosenbrock formula of order 2. See 
 Trapezoidal rule using a “free” interpolant. See 
 Implicit RungeKutta formula with a trapezoidal rule step as its
first stage and a backward differentiation formula of order 2 as its second
stage. See 
 Variablestep, variableorder (VSVO) solver using AdamsMoulton formulas, with the order varying between 1 and 12. Supports sensitivity analysis. See CVODE and CVODES. 
 Variablestep, variableorder (VSVO) solver using Backward Differentiation Formulas (BDFs) in fixedleading coefficient form, with order varying between 1 and 5. Supports sensitivity analysis. See CVODE and CVODES. 
 Variableorder, variablecoefficient solver using Backward Differentiation Formulas (BDFs) in fixedleading coefficient form, with order varying between 1 and 5. Supports sensitivity analysis. See IDA and IDAS. 
SelectedSolver
— Selected solver
name of a solver
This property is readonly.
Selected solver, returned as the name of a solver. If the value of
Solver
is "auto"
,
"stiff"
, or "nonstiff"
, then the
ode
object automatically chooses a solver and sets the value of
SelectedSolver
to the name of the chosen solver. If you
manually select a solver, then SelectedSolver
and
Solver
have the same value.
SolverOptions
— Solverspecific options
matlab.ode.options.*
object
Solverspecific options, specified as a matlab.ode.options.*
object. The ode
object automatically populates the
SolverOptions
property with an options object appropriate for
the selected solver. You can check available options for an existing
ode
object with the command
properties(F.SolverOptions)
. Specify options for the ODE problem
by changing property values of the matlab.ode.options.*
object. For
example, F.SolverOptions.OutputFcn = @odeplot
specifies an output
function that the solver calls after each successful time step.
When the Solver
property is "auto"
,
"stiff"
, or "nonstiff"
, the
SolverOptions
property is readonly.
This table summarizes the available options for each solver.
Options Object  Properties 



 




AbsoluteTolerance
— Absolute error tolerance
1e6
(default)  positive scalar  vector
Absolute error tolerance, specified as a positive scalar or vector. This tolerance
is a threshold below which the value of the solution becomes unimportant. If the
solution y
is less than AbsoluteTolerance
,
then the solver does not need to obtain any correct digits in y
.
For this reason, keep in mind the scale of the solution components when setting the
value of AbsoluteTolerance
.
If AbsoluteTolerance
is a vector, then it must be the same
length as the number of solution components. If AbsoluteTolerance
is a scalar, then the value applies to all solution components.
At each step, the ODE solver estimates the local error e
in the
i
th component of the solution. To be successful, the step must
have an acceptable error, as determined by both the relative and absolute error
tolerances:
e(i) <= max(RelativeTolerance*abs(y(i)),AbsoluteTolerance(i))
Data Types: single
 double
RelativeTolerance
— Relative error tolerance
1e3
(default)  positive scalar
Relative error tolerance, specified as a positive scalar. This tolerance measures
the error relative to the magnitude of each solution component. The relative error
tolerance controls the number of correct digits in all solution components, except
those smaller than AbsoluteTolerance
.
At each step, the ODE solver estimates the local error e
in the
i
th component of the solution. To be successful, the step must
have an acceptable error, as determined by both the relative and absolute error
tolerances:
e(i) <= max(RelativeTolerance*abs(y(i)),AbsoluteTolerance(i))
Data Types: single
 double
Object Functions
solve  Solve ODE over interval or at specified points 
solutionFcn  Construct function that interpolates ODE solution 
Examples
Solve Single ODE
Create an empty ode
object, and then specify values for the ODEFcn
and InitialValue
properties.
F = ode; F.ODEFcn = @(t,y) 2*t; F.InitialValue = 0;
Check which solver is selected for the problem, and then solve the equation over the time range [0 10]
.
F.SelectedSolver
ans = SolverID enumeration ode45
sol = solve(F,0,10)
sol = ODEResults with properties: Time: [0 0.2500 0.5000 0.7500 1 1.2500 1.5000 1.7500 2 2.2500 2.5000 2.7500 3 3.2500 3.5000 3.7500 4 4.2500 4.5000 4.7500 5 5.2500 5.5000 5.7500 6 6.2500 6.5000 6.7500 7 7.2500 7.5000 7.7500 8 8.2500 8.5000 8.7500 9 9.2500 9.5000 9.7500 10] Solution: [0 0.0625 0.2500 0.5625 1.0000 1.5625 2.2500 3.0625 4 5.0625 6.2500 7.5625 9 10.5625 12.2500 14.0625 16 18.0625 20.2500 22.5625 25 27.5625 30.2500 33.0625 36 39.0625 42.2500 45.5625 49 52.5625 56.2500 60.0625 64 68.0625 ... ] (1x41 double)
Plot the results.
plot(sol.Time,sol.Solution,"o")
Specify Jacobian Matrix for Stiff ODE
The Van der Pol oscillator equation is a secondorder differential equation. The equation includes a parameter $\mu $, and the equation becomes stiff when the value of $\mu $ is large.
$$\frac{{d}^{2}x}{d{t}^{2}}\mu (1{x}^{2})\frac{dx}{dt}+x=0$$
Using the substitutions ${\mathit{y}}_{1}=\mathit{x}$ and ${\mathit{y}}_{2}=\frac{\mathrm{dx}}{\mathrm{dt}}$ produces a system of two firstorder equations.
$$\begin{array}{l}\frac{d{y}_{1}}{dt}={y}_{2}\\ \frac{d{y}_{2}}{dt}=\mu (1{y}_{1}^{2}){y}_{2}{y}_{1}\end{array}$$
The Jacobian matrix for these equations is the matrix of partial derivatives of each equation with respect to both ${\mathit{y}}_{1}$ and ${\mathit{y}}_{2}$.
$$J=\left[\begin{array}{cc}\frac{\partial {f}_{1}}{\partial {y}_{1}}& \frac{\partial {f}_{1}}{\partial {y}_{2}}\\ \frac{\partial {f}_{2}}{\partial {y}_{1}}& \frac{\partial {f}_{2}}{\partial {y}_{2}}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 2\mu {y}_{1}{y}_{2}1& \mu (1{y}_{1}^{2})\end{array}\right]$$
Solve the Van der Pol oscillator using $\mu =1000$ and initial values of [2; 0]
by creating an ode
object to represent the problem.
Store the value of $\mu $ in the
Parameters
property.Specify the initial values in the
InitialValue
property.Specify the system of equations in the
ODEFcn
property, specifying three input arguments so that the value for $\mu $ is passed to the function.Specify a function that calculates the Jacobian matrix in the
Jacobian
property, specifying three input arguments so that the value for $\mu $ is passed to the function.
F = ode; F.Parameters = 1000; F.InitialValue = [2; 0]; F.ODEFcn = @(t,y,p) [y(2); p(1)*(1y(1)^2)*y(2)y(1)]; F.Jacobian = @(t,y,p) [0 1; 2*p(1)*y(1)*y(2)1 p(1)*(1y(1)^2)];
Display the ode
object. The SelectedSolver
property shows that the ode15s
solver was automatically chosen for this problem.
F
F = ode with properties: Problem definition ODEFcn: @(t,y,p)[y(2);p(1)*(1y(1)^2)*y(2)y(1)] Parameters: 1000 InitialTime: 0 InitialValue: [2x1 double] Jacobian: [1x1 odeJacobian] Solver properties AbsoluteTolerance: 1.0000e06 RelativeTolerance: 1.0000e03 Solver: auto SelectedSolver: ode15s Show all properties
Solve the system of equations over the time interval [0 3000]
by using the solve
method. Plot the first solution component.
S = solve(F,0,3000);
plot(S.Time,S.Solution(1,:),"o")
Solve System of DAEs with Mass Matrix
Consider this system of firstorder equations.
$$\begin{array}{l}{y}_{1}^{\prime}={y}_{1}{y}_{3}{y}_{2}\\ {y}_{2}^{\prime}={y}_{1}1\\ 0={y}_{1}+{y}_{2}+{y}_{3}\end{array}$$
The left side of the equations contain time derivatives, ${\mathit{y}}_{\mathit{n}}^{\prime}=\frac{{\mathrm{dy}}_{\mathit{n}}}{\mathrm{dt}}$. However, because the derivative for ${\mathit{y}}_{3}$ does not appear in the system, the equations define a system of differential algebraic equations. Rewriting the system in the form ${\mathrm{My}}^{\prime}=\mathit{f}\left(\mathit{t},\mathit{y}\right)$ shows a constant, singular mass matrix on the left side.
$$\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]\underset{}{\overset{\dot{}}{\left[\begin{array}{c}{y}_{1}^{\prime}\\ {y}_{2}^{\prime}\\ {y}_{3}^{\prime}\end{array}\right]}}=\left[\begin{array}{c}{y}_{1}{y}_{3}{y}_{2}\\ {y}_{1}1\\ {y}_{1}+{y}_{2}+{y}_{3}\end{array}\right]$$
Solve the system of equations using the initial values [1 1 2]
by creating an ode
object to represent the problem.
Specify the initial values in the
InitialValue
property.Specify the system of equations as an anonymous function in the
ODEFcn
property.Use an
odeMassMatrix
object to specify the constant, singular mass matrix in theMassMatrix
property.
F = ode;
F.InitialValue = [1 1 2];
F.ODEFcn = @(t,y) [y(1)*y(3)y(2);
y(1)1;
y(1)+y(2)+y(3)];
F.MassMatrix = odeMassMatrix(MassMatrix=[1 0 0; 0 1 0; 0 0 0],Singular="yes");
Display the ode object. The SelectedSolver
property shows that the ode15s
solver was automatically chosen for this problem.
F
F = ode with properties: Problem definition ODEFcn: @(t,y)[y(1)*y(3)y(2);y(1)1;y(1)+y(2)+y(3)] InitialTime: 0 InitialValue: [1 1 2] Jacobian: [] MassMatrix: [1x1 odeMassMatrix] Solver properties AbsoluteTolerance: 1.0000e06 RelativeTolerance: 1.0000e03 Solver: auto SelectedSolver: ode15s Show all properties
Solve the system of equations over the time interval [0 10]
by using the solve
method. Plot all three solution components.
S = solve(F,0,10); plot(S.Time,S.Solution,"o") legend("y_1","y_2","y_3",Location="southeast")
Examine Parameter Sensitivity
Solve an ODE system with two equations and two parameters, and perform sensitivity analysis on the parameters.
Create an ode
object to represent this system of equations.
$$\begin{array}{l}\frac{{\mathrm{dy}}_{1}}{\mathrm{dt}}={\mathit{p}}_{1}{\mathit{y}}_{1}{\mathit{y}}_{2}\\ \frac{{\mathrm{dy}}_{2}}{\mathrm{dt}}={\mathit{p}}_{2}{\mathit{y}}_{2}\end{array}$$
Specify the initial conditions as ${\mathit{y}}_{1}\left(0\right)=2$ and ${\mathit{y}}_{2}\left(0\right)=3$, and parameter values of ${\mathit{p}}_{1}=0.05$ and ${\mathit{p}}_{2}=1.5$. To enable sensitivity analysis of the parameters, set the Sensitivity
property of the ode
object to an odeSensitivity
object.
p = [0.05 1.5]; F = ode(ODEFcn=@(t,y,p) [p(1)*y(1)y(2); p(2)*y(2)], ... InitialValue=[2 3], ... Parameters=p, ... Sensitivity=odeSensitivity)
F = ode with properties: Problem definition ODEFcn: @(t,y,p)[p(1)*y(1)y(2);p(2)*y(2)] Parameters: [0.0500 1.5000] InitialTime: 0 InitialValue: [2 3] Sensitivity: [1x1 odeSensitivity] Solver properties AbsoluteTolerance: 1.0000e06 RelativeTolerance: 1.0000e03 Solver: auto SelectedSolver: cvodesNonstiff Show all properties
Because the equations are nonstiff and sensitivity analysis is enabled, the ode
object automatically chooses the cvodesNonstiff
solver for this problem.
Solve the ODE over the time interval [0 5]
, and plot the solution for each component.
S = solve(F,0,5)
S = ODEResults with properties: Time: [0 2.9540e09 2.9543e05 2.2466e04 4.1978e04 0.0024 0.0139 0.0255 0.0484 0.0713 0.1363 0.2014 0.2992 0.3970 0.4948 0.6661 0.8373 1.0085 1.1798 1.3510 1.5222 1.8047 2.0872 2.3697 2.6522 2.9347 3.2172 3.4997 3.7822 ... ] (1x34 double) Solution: [2x34 double] Sensitivity: [2x2x34 double]
plot(S.Time,S.Solution(1,:),"o",S.Time,S.Solution(2,:),"o") legend("y1","y2")
The values in S.Sensitivity
are partial derivatives of the equations with respect to the parameters. To examine the effects of the parameter values during the integration, plot the sensitivity values.
figure hold on plot(S.Time,squeeze(S.Sensitivity(1,1,:)),"o") plot(S.Time,squeeze(S.Sensitivity(1,2,:)),"o") plot(S.Time,squeeze(S.Sensitivity(2,1,:)),"o") plot(S.Time,squeeze(S.Sensitivity(2,2,:)),"o") legend("p1,eq1","p2,eq1","p1,eq2","p2,eq2") hold off
Use Callback Function to Restart Integration of Bouncing Ball
Consider a ball thrown upward with an initial velocity $\mathrm{dy}/\mathrm{dt}$. The ball is subject to acceleration due to gravity aimed downward, so its acceleration is
$$\frac{{d}^{2}y}{d{t}^{2}}=g.$$
Rewriting the equation as a firstorder system of equations with the substitutions ${\mathit{y}}_{1}=\mathit{y}$ and ${\mathit{y}}_{2}=\frac{\mathrm{dy}}{\mathrm{dt}}$ yields
$$\begin{array}{l}{y}_{1}^{\prime}={y}_{2}\\ {y}_{2}^{\prime}=g.\end{array}$$
Solve the equations for the position ${\mathit{y}}_{1}$ and velocity ${\mathit{y}}_{2}$ of the ball over time.
Define Equations and Initial Conditions
Create a function handle for the firstorder system of equations that accepts two inputs for (t,y)
. Use the value $\mathit{g}=9.8$ for the acceleration due to gravity.
dydt = @(t,y) [y(2); 9.8];
Next, create a vector with the initial conditions. The ball starts at position ${\mathit{y}}_{1}=3$ at $\mathit{t}=0$ as it is thrown upward with initial velocity ${\mathit{y}}_{2}=20$.
y0 = [3 20];
Model Ball Bounces as Events
The ball initially travels upward until the force due to gravity causes it to change direction and head back down to the ground. If you solve the equations without more consideration, then the ball falls back downward forever without striking the ground. Instead, you can use an event function to detect when the position of the ball goes to zero where the ground is located. Because the solution component ${\mathit{y}}_{1}=\mathit{y}$ is the position of the ball, the event function tracks the value of ${\mathit{y}}_{1}$ so that an event triggers whenever ${\mathit{y}}_{1}=0$.
Create a function handle for the event function that accepts two inputs for (t,y)
.
bounceEvent = @(t,y) y(1);
When the ball strikes the ground, its direction changes again as it heads back upwards with a new (smaller) initial velocity. To model this situation, use a callback function along with the event function. When an event triggers, the ODE solver invokes the callback function. The callback function resets the position and initial velocity of the ball so that the integration can restart with the correct initial conditions. When an event occurs, the callback function sets the position ${\mathit{y}}_{1}=0$ and attenuates the velocity by a factor of $0.9$ while reversing its direction back upward. Define a callback function that performs these actions.
function [stop,y] = bounceResponse(t,y) stop = false; y(1) = 0; y(2) = 0.9*y(2); end
(The callback function is included as a local function at the end of the example.)
Create an odeEvent
object to represent the bouncing ball events. Specify Direction="descending"
so that only events where the position is decreasing are detected. Also, specify Response="callback"
so that the solver invokes the callback function when an event occurs.
E = odeEvent(EventFcn=bounceEvent, ... Direction="descending", ... Response="callback", ... CallbackFcn=@bounceResponse)
E = odeEvent with properties: EventFcn: @(t,y)y(1) Direction: descending Response: callback CallbackFcn: @bounceResponse
Solve Equations
Create an ode
object for the problem, specifying the equations dydt
, initial conditions y0
, and events E
as property values.
F = ode(ODEFcn=dydt,InitialValue=y0,EventDefinition=E);
Integrate the equations over the time interval [0 30]
by using the solve
method. Specify Refine=8
to generate 8 points per step. The resulting object has properties for the time and solution, and because events are being tracked, the object also displays properties related to the events that triggered during the integration.
S = solve(F,0,30,Refine=8)
S = ODEResults with properties: Time: [0 0.0038 0.0075 0.0113 0.0151 0.0188 0.0226 0.0264 0.0301 0.0490 0.0678 0.0867 0.1055 0.1243 0.1432 0.1620 0.1809 0.2751 0.3692 0.4634 0.5576 0.6518 0.7460 0.8402 0.9344 1.3094 1.6844 2.0594 2.4344 2.8094 3.1844 ... ] (1x468 double) Solution: [2x468 double] EventTime: [4.2265 8.1607 11.7015 14.8882 17.7563 20.3375 22.6606 24.7514 26.6331 28.3267 29.8509] EventSolution: [2x11 double] EventIndex: [1 1 1 1 1 1 1 1 1 1 1]
Plot Results
Plot the position ${\mathit{y}}_{1}$ of the ball over time, marking the initial position with a green circle and events with red circles.
plot(S.Time,S.Solution(1,:),"") hold on plot(S.EventTime,S.EventSolution(1,:),"ro") plot(0,y0(1),"go") hold off ylim([0 25]) xlabel("Time") ylabel("Position y_1")
Local Functions
function [stop,y] = bounceResponse(t,y) stop = false; y(1) = 0; y(2) = 0.9*y(2); end
Version History
Introduced in R2023b
MATLAB Command
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
 América Latina (Español)
 Canada (English)
 United States (English)
Europe
 Belgium (English)
 Denmark (English)
 Deutschland (Deutsch)
 España (Español)
 Finland (English)
 France (Français)
 Ireland (English)
 Italia (Italiano)
 Luxembourg (English)
 Netherlands (English)
 Norway (English)
 Österreich (Deutsch)
 Portugal (English)
 Sweden (English)
 Switzerland
 United Kingdom (English)