Download Using ODE45-Matlab-Mannual and more Slides Matlab skills in PDF only on Docsity!
ENGR210 Using ODE45 1
Bucknell University Using ODE45 MATLAB Help
MATLAB's standard solver for ordinary differential equations (ODEs) is the function ode45. This function implements a Runge-Kutta method with a variable time step for efficient computation. ode45 id designed to handle the following general problem
d dt
t t (^) o o
y = f ( , y ) y ( ) = y [1]
where t is the independent variable (time, position, volume) and y is a vector of dependent variables (temperature, position, concentrations) to be found. The mathematical problem is specified when the vector of functions on the right-hand side of Eq. [1], f ( , t y ), is set and the initial conditions, y =
y o at time t o , are specified.
The notes here apply to versions of MATLAB above 5.0 and cover the basics of using the function ode45. For more information on this and other ODE solvers in MATLAB, see the on-line help.
Contents:
Syntax for ode45 ....................................................................................................................... 2 Integrating a single, first-order equation...................................................................................... 3 Getting the solution at particular values of the independent variable........................................... 4 Integrating a set of coupled first-order equations......................................................................... 4 Integrating a second-order initial-value problem (IVP)................................................................ 7 Integrating an N th-order initial-value problem............................................................................ 8 Changing model parameters........................................................................................................ 9 Integrating a second-order boundary-value problem (BVP)........................................................ 11 Setting options in ode45 .......................................................................................................... 12 Going beyond ode45 ................................................................................................................. 13
ENGR210 Using ODE45 2
Syntax for ode
ode45 may be invoked from the command line via
[t,y] = ode45('fname', tspan, y0, opts)
where
fname name of the function Mfile used to evaluate the right-hand-side function in Eq. [1] at a given value of the independent variable and dependent variable(s) (string). The function definition line usually has the form
function dydt = fname(t,y)
The output variable ( dydt ) must be a vector with the same size as y. Note that the independent variable ( t here) must be included in the input argument list even if it does not explicitly appear in the expressions used to generate dydt.
tspan 2-element vector defining the range of integration ( [to tf] ) though variations are possible.
y0 vector of initial conditions for the dependent variable. There should be as many initial conditions as there are dependent variables.
opts a MATLAB structure variable that allows you to control the details of computation (if you want to). This argument is optional and, if not provided, ode45 will use default values (see the examples below).
t Value of the independent variable at which the solution array ( y ) is calculated. Note that by default this will not be a uniformly distributed set of values.
y Values of the solution to the problem (array). Each column of y is a different dependent variable. The size of the array is length(t) -by- length(y0)
Specific examples of using ode45 now follow. Mfiles for these examples are in the body of this document and should also be available in the folder that contains this document.
ENGR210 Using ODE45 4
Getting the solution at particular values of the independent variable
ode45 uses a variable-step-length algorithm to find the solution for a given ODE. Thus, ode varies the size of the step of the independent variable in order to meet the accuracy you specify at any particular point along the solution. If ode45 can take "big" steps and still meet this accuracy, it will do so and will therefore move quickly through regions where the solution does not "change" greatly. In regions where the solution changes more rapidly, ode45 will take "smaller" steps. While this strategy is good from an efficiency or speed point of view, it means that the solution does not appear at a fixed set of values for the independent variable (as a fixed-step method would) and sometimes the solution curves look a little ragged.
The simplest way to improve on the density of solution points is to modify the input tspan from a 2-element vector to an N -element vector via something like
>> tspan = linspace(to,tf,500)’;
and use this new version in the input list to ode.
Smoother curves can also be generated by interpolation (spline interpolation usually works nicely). For example, if you wanted a smoother result from the solution for the tank-fill problem, you might do the following
>> ti = linspace(tspan(1),tspan(2),300); (300 points - you could use more) **>> hi = spline(t,h,ti);
plot(t,h,’o’,ti,hi);**
The interpolated curve smoothes out the rough edges caused by simply connecting the data points (which is what plot does) and so makes the graph more appealing, in a visual sense.
Integrating a set of coupled first-order equations
Chemical-kinetics problems often lead to sets of coupled, first-order ODEs. For example, consider the reaction network
A ↔ B → C [4]
Assuming a first-order reaction-rate expression for each transformation, material balances for each species lead to the following set of ODEs:
dA dt
k A k B
dB dt
k A k B k B
dC dt
k B
1 2
1 2 3
3
[5]
with the initial conditions, (^) A ( ) 0 = A (^) o , B ( ) 0 = B (^) o , C ( ) 0 = Co. Since the equations are coupled, you
cannot solve each one separately and so must solve them simultaneously.
The system in Eq. [5] can be put in the standard form for ode45 (Eq. [1]) by defining the vectors y , y o and f as
ENGR210 Using ODE45 5
y = y y f y
A
B
C
A
B
C
t
k y k y k y k k y k y
o
o o o
1 1 2 2 1 1 2 3 2 3 2
[6]
Solving the system represented by Eq. [6] is a simple extension of what was done for solving a single equation. We'll demonstrate the solution for the following situation
k 1 (^) = 5 k (^) 2 = 2 k (^) 3 = 1 A (^) o = 1 Bo = Co = 0
Step 1: Write a function Mfile to evaluate the right-hand-side expression
The primary difference here, compared to the single-equation case, is that the input variable y will be a vector. The first element of y represents the concentration of species A at a time t , and the second and third elements representing the concentrations of species B and C, respectively, at the same time, t. This ordering of variables is defined by Eq. [6]. There is no "right" order to the variables but whatever order you do choose, use it consistently. We'll call the Mfile react.m. It looks like this:
function dydt = react(t,y) % Solve the kinetics example
dydt = zeros(size(y));
% Parameters - reaction-rate constants
k1 = 5; k2 = 2; k3 = 1;
A = y(1); We'll be explicit about it here though you can do B = y(2); the calculations directly with the y -values. C = y(3);
% Evaluate the RHS expression
dydt(1) = -k1A + k2B; dydt(2) = k1A - (k2+k3)B; dydt(3) = k3B;*
% eof - react.m
Note that the input arguments must be t and y (in that order) even though t is not explicitly used in the function.
Step 2: Use ode45 to solve the problem
No time interval is given so we'll pick one (0 to 4) and see what the solution looks like. If a longer or shorter interval is needed, we can simply re-execute the function with a new value for the ending time. Following the outline for the single-equation problem, the call to ode45 is,
>> [t,y] = ode45('react',[0 4],[1 0 0]);
Note that the initial condition is provided directly in the call to ode45. You could also have defined a variable y0 prior to the call to ode45 and used that variable as an input.
ENGR210 Using ODE45 7
Integrating a second-order initial-value problem (IVP)
A mass-spring-dashpot system can be modeled via the following second-order ODE
y^ ˙˙^ + cy ˙^ + ω^2 y = g t ( ) y ( ) 0 = y (^) o , v ( ) 0 = y ˙ ( ) 0 = vo [7]
In this model, c represents a retarding force proportional to the velocity of the mass, ω^2 is the natural frequency of the system and g ( t ) is the forcing (or input) function. The initial conditions are the initial position ( y o ) and initial velocity ( v o ).
ode45 is set up to handle only first-order equations and so a method is needed to convert this second- order equation into one (or more) first-order equations which are equivalent. The conversion is accomplished through a technique called "reduction of order". We'll illustrate the solution for the particular set of conditions
c = 5 ω = 2 y ( ) 0 = 1 v ( ) 0 = 0 g t ( ) =sin( ) t
Step1: Define the components of a vector p = [ p p 1 2 ]T^ as follows:
p y p y
1 2
= ˙^
[8]
Step 2: Form the first derivatives of each of the components of p
Using the given differential equation, we can write a system of first-order equations as
˙ ˙ ˙ ˙˙ (^) ( ) ˙
( )
p y p p y g t cy y g t cp p
1 2 2
2
2
2 1
ω ω
[9]
In writing the expression for the second component, we've used the governing ODE (Eq. [7]).
Step 3: Cast the problem in the format needed to use ode.
p ( , )
p p p
= = f p
d dt
d dt
p p
p g t cp p
f t f t
(^1) t 2
2 2
2 1
1 ω 2
[10]
Step 4: Collect the initial conditions into a single vector
p ( ) p
1 2
o
o o
p p
y y
y v
[11]
Step 5: Apply ode45 to solve the system of equations
The Mfile for the RHS function for this problem will be called spring.m. Here it is:
function pdot = spring(t,p) % Spring example problem
ENGR210 Using ODE45 8
% Parameters - damping coefficient and natural frequency
c = 5; w = 2;
g = sin(t); % forcing function pdot = zeros(size(p)); pdot(1) = p(2); pdot(2) = g - cp(2) - (w^2)p(1);**
% eof - spring.m
The call to ode45 is, for a solution interval of 0 to 20,
>> p0 = [1 0]; (initial position and velocity) >> [t,p] = ode45('spring',[0 20],p0);
Step 6: Look at the results
If you wanted to look at only the displacement, you'd want to look at the first column of p (see the definition of p in the first step, Eq. [8]). Hence, you would give the command
>> plot(t,p(:,1))
An interesting plot for these sorts of problems is the phase-plane plot, a plot of the velocity of the mass versus its position. This plot is easily created from your solution via
>> plot(p(:,1),p(:,2))
Phase-plane plots are useful in analyzing general features of dynamic systems.
Integrating an N th-order initial-value problem
To use ode45 to integrate an Nth-order ODE, you simply continue the process outlined in the section on integrating a 2nd-order ODE. The first element of the vector p is set to the dependent variable and then subsequent elements are defined as the derivatives of the dependent variable up to one less than the order of the equation. Finally, the initial conditions are collected into one vector to give the format presented in Eq. [1].
For example, the 4th-order equation
a
d y dx
b
d y dx
c
d y dx
d
dy dx
ey
4 4
3 3
2
would generate the first-order system
d dt
p p p p
p p p bp cp dp ep a
1 2 3 4
2 3 4 4 3 2 1
[13]
which, along with an appropriate set of initial conditions would complete the set-up for ode.
ENGR210 Using ODE45 10
if isempty(C_SPRING) % how to get out break end
W_SPRING = input('Natural frequency [w]: ');
[t,p] = ode45('spring2',[0 20],[1 0],[],C_SPRING,W_SPRING);
plot(p(:,1),p(:,2)) title('Phase-plane plot of results') xlabel('position') ylabel('velocity')
end
% eos - dospring.m
Note the additions to the call to ode45. First, a placeholder for the “options” input is inserted (an empty array) so that “default” options are used. Then, the parameters for the model are provided in the order that they appear in the definition line of the RHS function.
Try the script out and modify it (e.g., you could add the frequency and/or amplitude of the forcing function as something to be changed).
ENGR210 Using ODE45 11
Integrating a second-order boundary-value problem (BVP)
ode45 was written to solve initial-value problems (IVPs). Hence it cannot be (directly) used to solve, for example, the following problem derived from a model of heat-transfer in a rod:
d
2
y
dx^2
− y = 0 y (0) = 1
dy
dx x = 1
= 0 [14]
since the value of the derivative at x = 0 is not specified (it is known at x = 1, though). Equation [14] is a boundary-value problem (BVP) and is common in models based on transport phenomena (heat transfer, mass transfer and fluid mechanics).
All is not lost because one way to solve a BVP is to pretend it is an IVP. To make up for the lack of knowledge of the derivative at the initial point, you can guess a value, do the integration and then check yourself by seeing how close you are to meeting the conditions at the other end of the interval. When you have guessed the right starting values, you have the solution to the problem. This approach is sometimes called the "shooting method" by analogy to the ballistics problem of landing an artillery shell on a target by specifying only it's set-up (powder charge and angle of the barrel).
Step 1: Set up the problem so that ode45 can solve it
Using the approach of turning a second order equation into a pair of coupled first-order equations, we have
d dx
p p
p p v
1 2
2 1
p ( ) [15]
where v has been used to represent the (unknown) value of the derivative at x = 0. The Mfile used to evaluate the RHS is as follows
function dpdx = hotrod(x,p) % Hot-rod problem illustrating the shooting method
dpdx = zeros(size(p));
dpdx(1) = p(2); dpdx(2) = p(1);
% eof - hotrod.m
Step 2: Guess a value of the initial slope and integrate to x = 1
The problem will be iterative so it's not likely that the first guess will be right. From the physics of the problem, the end of the rod (at x = 1) will be colder than the place we are starting from ( x = 0) and so we'll guess a negative value for the initial slope.
**>> v = -1;
[x,p] = ode45('hotrod',[0 1],[1 v]);**
The value of the derivative at x = 1 is the last value in the second column of p (why?). Thus, we can check the accuracy of the first guess via
ENGR210 Using ODE45 13
Going beyond ode
The solver ode45 is not the be-all and end-all of ODE-solvers. While ode45 should be your first choice for integration, there are problems that the function performs poorly on or even fails on. In such cases, there are fallback solvers that are available. All these solvers use the same syntax as ode45 (see page 2) but have options for handling more difficult or sophisticated problems.
Here are some suggestions for handling non-standard ODE problems:
- If accuracy you desire is not obtainable via ode45 , try the function ode113. This solver uses a variable order method that may be able to improve over what ode45 does.
- If ode45 is taking too long to compute a solution, your problem may be “stiff” (i.e., it involves a system with a wide range of time constants). Try the function ode15s.
- If your system of equations has the form
M
y f y
d dt
= ( , t )
where M is a (typically non-singular) matrix, try the function ode15s.
You’ll find more information on these function in the on-line help and documentation. For example, try the on-line function reference (available through the command helpdesk ) on any of the solvers noted above.
