11 Assigned: November 2011
Fall 2011
Due: 18 November 2011
E7 Laboratory Assignment 12
This lab covers the simulation of ordinary dierential equations (ODEs) and statistics. The range of
ODE simulations provided in this lab should prepare you to simulate any ODE in the future. The
statistics problems provided should equip you with valuable techniques to analze data with.
Note: You will use MATLABs publish command (accessible in the Editor window via File ->
Publish Filename) to generate a printable document that contains your code and your answers.
1. You must turn in a printed copy to the homework boxes in room 1109 Etcheverry. Further, you will
need to use the command type to display your functions in your main script. You will not get
credit for the functions you have written unless you display them! Be sure to turn your printout
into the box corresponding to the lab section in which you are registered. The due date is 12:00 pm
(noon) on November 18, 2011.
2. Your primary m-le, which is used to generate the HTML document, and all of the function m-les
that you create must also be uploaded to bSpace (under the Assignments tab) by noon on November
18, 2011.
Please (i) have your full name, student ID number, and lab section number at the head of the le and
(ii) name this le lastname middlename firstname lab12.m.
Use the clear; format compact commands at the beginning of your script to clear the variables from
the base workspace and format appropriately to reduce the length of the output. If your comments
exceed the page border, highlight the comment and go to Text->Wrap Comment.
In several problems, you are asked to execute a number of given lines, and explain (in a comment before
the line) what the command does. Some lines might produce an error (this is deliberate, on our part).
In that case, after you understand what is the issue, comment out that line (using the % symbol) in
your script le and explain why the error occurred.
Part A: Introduction to ODEs
1. Warm-up problem: using ode45 to solve 1st order ODEs.
Start by creating anonymous functions for the ODEs listed below. Use a time span of 0 to 10
seconds, i.e., tspan = [0 10]. Note that by default ode45 will use an adaptive step size when
simulating the ODE. Use initial conditions of 1 for each element in the dependent variable vector
x and refer to the lecture notes for the use of ode45.
Plot the dependent variable vector x as a function of the independent variable (here, time t), i.e.,
plot(t,x). If an analytical solution is provided, include it as well in the plot, and be sure to label
each solution curve within the legend.
In the following three cases, the dependent variable is a scalar:
1 of 5
Assigned: 11 November 2011
Fall 2011
Due: 18 November 2011
x = cos t (analytical solution: x = 1 + sin t)
x = x (analytical solution x = et )
x = x (analytical solution x = et )
Another option is for the function handle that represents the ODE to contain other inputs besides
the independent and dependent variables. To implement this, create a function handle for the
ODE given below with an additional input argument pVal that represents the value of the scalar
variable p.
x = px (analytical solution x = ept )
Now, use ode45 to simulate this ODE for the following three cases: p = 0.1, p = -0.5, and p =
-0.1. In order to do this, you need to create a function handle, called @fp with arguments t and
x, that will call another function handle, called odefun, that denes the created with right-hand
side of the ODE. An example of this for p = 0.1 is shown below:
begin code
1
2
odefun = @(t,x,p) p*x;
fp = @(t,x) odefun(t,x, 0.1);
end code
Simulate all three cases and plot the state vectors of each in the same gure with a legend.
In the next case, the dependent variable is a 2 1 vector. Simulate this ODE and plot both of
the state vectors in the same gure with a legend.
x=
x2
x1
=
x1
x2
2. Warm-up problem: usage of ode45 to solve higher-order ODEs.
Start by creating anonymous functions for the ODEs listed below. You will have to transform the
higher-order (greater than 1) ODEs to a series of 1st order ODEs. The order of an ODE is dened
by the order of the highest order derivative contained in the ODE. Use a time span of 0 to 10
seconds. Use initial conditions of 1 for each element in the state vector. Plot the state vector x
as a function of time t, i.e., plot(t,x).
This is a 2nd order system.
x+x+x=0
This is a forced 2nd order system.
x + x + x = f (t) where f (t) = sin(t)
You should be able to use the techniques used in this problem to solve ODEs of any order.
2 of 5
Assigned: November 11 2011
Fall 2011
Due: 18 November 2011
Part B: Simulating Population Dynamics
3. In this problem, you will use ode45 to simulate a population dynamics problem, sometimes referred to as a predator-pray (e.g., fox-rabbit) problem. These problems predict the growth (or
shrinkage) in the population of two species (a predator and a prey species) based on some
assumed model of their interaction.
Let the sizes of the prey and predator populations be denoted x and y , respectively, and assume
that the population growth of the two species is described by the rst-order ODEs:
Prey population growth, x = x(a1 a2 y ) preyyear1
Predator population growth, y = y (a3 a4 x) predatorsyear1
The parameters entering the preceding equations take the following values:
Prey growth rates, a1 = 0.55 year1 , a2 = 0.034 predators1 year1
Predator growth rates, a3 = 0.9 year1 , a4 = 0.011 prey1 year1
Clearly, the two equations are coupled, therefore need to be solved as a system (as opposed to be
solved separately).
Start by creating anonymous functions for the system of ODEs. Then, solve the system using a
time span of 0 to 100 years and initial conditions of 100 prey and 5 predators. Finally, generate a
single plot of the prey population variable x and the predator population variable y as a function
of time t. The version of the expected plot is included below for reference.
a1 = 0.550, a2 = 0.034, a3 = 0.900, a4=0.011
180
Prey
Predators
160
140
Population
120
100
80
60
40
20
0
0
10
20
30
40
50
Time (Years)
60
70
80
Figure 1: Predator-Prey Populations
3 of 5
90
100
Assigned: 11 November 2011
Fall 2011
Due: 18 November 2011
Part C: Simulating Physical Systems
4. In this problem, you will use ode45 to simulate the motion of the double-pendulum illustrated in
the gure below.
g
L1
L2
y
x
Figure 2: Double Pendulum
The equations of motion of the double pendulum are as follows:
=
m2 L1 2 sin( ) cos( )+gm2 sin( ) cos( )m2 L2 2 sin( )(m1 +m2 )g sin()
L1 (m1 +m2 )m2 L1 cos( )2
=
m2 L2 2 sin( ) cos( )+g sin() cos( )(m1 +m2 )+L1 2 sin( )(m1 +m2 )g sin( )(m1 +m2 )
L2 (m1 +m2 )m2 L2 cos( )2
Here and are the angles of the two rods of the pendulum from the vertical conguration, as in
the gure. Also, m1 , m2 are the two masses (located at the ends of the two rods of length L1 and
L2 , respectively), and g is the gravitational acceleration. Assume that the preceding parameters
take the following values:
Gravitational acceleration, g = 9.81 ms2
Length of rod 1, L1 = 1 m
Length of rod 2, L2 = 1 m
Mass of joint 1, m1 = 1 kg
Mass of joint 2, m2 = 1 kg
Two cases of the double pendulum problems will be analyzed, with initial conditions dened as
follows:
4 of 5
Assigned: 11 November 2011
Fall 2011
Due: 18 November 2011
(a) Case 1
Angular
Angular
Angular
Angular
position of rod 1, = 0.5 rad
velocity of rod 1, = 0 rads1
position of rod 2, = 0.5 rad
velocity of rod 2, = 0 rads1
(b) Case 2
Angular
Angular
Angular
Angular
position of rod 1, = 0.001 rad
velocity of rod 1, = 0 rads1
position of rod 2, = + 0.001 rad
velocity of rod 2, = 0 rads1
Start by creating an anonymous function, called dpend, for the double pendulum ODEs that you
will subsequently use to call ode45.
1
Use a time span of 0 to 20 seconds, but specify that the outputs of ode45 should be every 30
seconds. To do this dene the time span as an evenly spaced vector from 0 to 20 with increments
1
of 30 , i.e., dene tspan = 0:1/30:20. Note when solving the problem, ode45 will still use an
adaptive step size to control error, but will output values at the intervals specied by the time
span vector.
Below is a script to animate the double pendulum. Include this in your script le so that the nal
gure is displayed in your published le. Note that handle graphics would be a better way to do
this animation. In order for the following script to run properly, the outputs of your call to ode45
should be the time vector, t, and the state vector, x. The state vector should have the form:
x = [ ].
begin code
1
2
3
4
5
6
7
8
9
10
11
12
13
% Position of mass 1 with respect to the origin
x1 = L1*sin(x(:,1)); y1 = -L1*cos(x(:,1));
% Position of mass 2 with respect to the origin
x2 = x1+L2*sin(x(:,3)); y2 = y1-L2*cos(x(:,3));
for i = 1:length(t)
clf; plot([0 x1(i)],[0 y1(i)],b); hold on;
plot(x1(i),y1(i),r.);
plot([x1(i) x2(i)],[y1(i) y2(i)],b);
plot(x2(i),y2(i),r.);
plot(x2(1:i),y2(1:i),k--);
title(sprintf(Time: %.1f seconds,t(i))); xlabel(x); ylabel(y);
axis equal; axis([-2,2,-2,2]); hold off; pause(1/30); drawnow;
end
end code
5 of 5