E7 Fall 2019, UC Berkeley Homework Assignment 12 Ordinary Differential Equations The purpose of this lab is to introduce you to ODEs and ODE solvers in MATLAB. Note: Use the templates provided in the...

E7 Fall 2019, UC Berkeley
Homework Assignment 12
Ordinary Differential
Equations
The purpose of this lab is to introduce you to ODEs and ODE solvers in MATLAB.
Note: Use the templates provided in the assignment download to complete this assignment.
The template has a very specific format to enable auto-grading, and if the format is altered
your assignment may not be properly graded (this will be apparent to you).
The due time is 11:59pm (midnight) December 11, 2019. NO LATE HOME-
WORK WILL BE ACCEPTED. Please upload the following files through the bCourses
website:
• HW12 ODE.m
• HW12 ODE.pdf (created using the publish command)
• All function files you are asked to create or are necessary to run your code.
• All generated plot files specified in the questions.
• Generated movie file specified.
Directions to upload files can be found here.
NOTE: This homework is designed in a sequential manner and gradually builds some con-
cepts as the problems progress. Some latter problems might have dependencies on initial
problems.
1. The Lorenz attractor, named after Edward N. Lorenz, was derived as a simplified model
of convection cu
ents in the atmosphere. It is known for having chaotic, unpredictable
ehavior. The Lorenz attractor can be described with the following three coupled
ordinary differential equations:
dx
dt
= σ(y − x)
dy
dt
= x(ρ− z) − y
dz
dt
= xy − βz
1
https:
guides.instructure.com/m/4212/l/54353-how-do-i-upload-a-file-to-my-assignment-submission
E7 Fall 2019, UC Berkeley
where t is the time, x, y, and z are the state variables, and σ, β, and ρ are scala
parameters. σ is called the Prandtl number and ρ is called the Rayleigh number.
Usually σ = 10, β = 8/3 and ρ is varied. The system exhibits chaotic behavior fo
ρ = 28, but displays knotted periodic o
its for other values of ρ.
The ODEs that describe the Lorentz attractor can be written in vector form Ẏ = f(Y )
d
dt
 Y1Y2
Y3
 =
 σ (Y2 − Y1)Y1 (ρ− Y3) − Y2
Y1 Y2 − β Y3
 (1)
Here Y1 is x, Y2 is y and Y3 is z. Write a function named LorenzODE with the heade
line:
function dotY = LorenzODE(Y, sigma, beta, rho)
whose first input argument Y is a 3x1 column vector and sigma, beta, and rho are
constants, and returns the right hand side of Eq. (1) (also a 3x1 column vector).
Use ode45 to numerically integrate the ODE Eq. (1) for: σ = 10, β = 8/3, and ρ = 14,
tspan = [0,20], and the following set of initial conditions:
(a) Y0 = [0;1;1]
(b) Y0 = [1;1;1]
(c) Y0 = [1;-1;1]
(d) Y0 = [1;-1;-1]
(e) Y0 = [1;-1;-10]
Remember that first input argument to ode45 should be an anonymous function handle,
such as
Lorenz = @(t,Y) LorenzODE(Y, sigma, beta, rho)
You can use the commands for plotting:
p lo t3 (Y( : , 1 ) ,Y( : , 2 ) ,Y( : , 3 ) ) ;
view ( 0 , 0 ) ;
hold on ;
to do 3D plots of the solution (Y1(t), Y2(t), Y3(t)) for all five initial conditions on the
same figure (using different values of color). Add the initial conditions as legend entries.
Save the figure using print command as ’Lorenz.png’.
2. In this problem, you will use ode45 to simulate the motion of the double pendulum
illustrated in the figure below.
The equations of motion of the double pendulum are as follows:
2
E7 Fall 2019, UC Berkeley
Figure 1: Double Pendulum System
• θ̈ = −m2L1θ̇
2 sin(θ−β) cos(θ−β)+gm2 sin(β) cos(θ−β)−m2L2β̇2 sin(θ−β)−(m1+m2)g sin(θ)
L1(m1+m2)−m2L1 cos2(θ−β)
• β̈ = m2L2β̇
2 sin(θ−β) cos(θ−β)+g sin(θ) cos(θ−β)(m1+m2)+L1θ̇2 sin(θ−β)(m1+m2)−g sin(β)(m1+m2)
L2(m1+m2)−m2L2 cos2(θ−β)
Here θ and β are the angles of the two rods of the pendulum from the vertical con-
figuration, as in the figure. 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 m·s−2
• 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
Use initial conditions defined as follows:
• Angular position of rod 1, θ = 0.5 rad
• Angular velocity of rod 1, θ̇ = 0 rad·s−1
• Angular position of rod 2, β = 0.5 rad
• Angular velocity of rod 2, β̇ = 0 rad·s−1
Notice that there are two equations provided but have double derivatives (angula
acceleration). We can, however, simplify them by identifying that θ̇ and β̇ can be set
as two variables and write two equations for them, i.e.
d
dt
θ = θ̇
3
E7 Fall 2019, UC Berkeley
d
dt
β = β̇
Thus the whole system can be written as Ẏ = f(Y ) as in the previous question.
Start by creating an anonymous function, called dpend, for the double pendulum ODEs
that you will subsequently use to call ode45. The function returns the value of f(Y) fo
an input of Y. Here Y = [θ; θ̇; β; β̇].
Use a time span of 0 to 20 seconds, but specify that the outputs of ode45 should be
every 1
30
seconds. To do this define the time span as an evenly spaced vector from 0
to 20 with increments of 1
30
, i.e., define tspan = 0:1/30:20. Note when solving the
problem, ode45 will still use an adaptive step size to control e
or, but will output
values at the intervals specified by the time span vector.
Define a function Y=DoublePendulum(I,t,C). Here I are the initial condition vecto
I = [θinitial; θ̇initial; βinitial; β̇initial], t is the time span vector tspan, C is the constant
vector C = [L1;L2;m1;m2; g] and Y is the output vector with each column Yi defined
as Yi = [θi; θ̇i; βi; β̇i] at time step ti. The total number of column vectors will be equal
to the t vector, i.e. every column describes the angular position and angular velocities
of the system.
Remembering that θ and β are angular positions, use clf and drawnow commands to
create an animation and save that animation as ’DoublePendulum.avi’ movie. Remem-
er to use suitable axis limits so the limits do not change from frame to frame. You
have already done this for a previous homework. You can display the pendulums as ’o’
and rods with with dashed lines.
3. In this problem, you will come up with your own integration functions using Euler’s
method and Runge-Kutta 4 (RK4) methods. Consider the following system of equa-
tions:
̈ − rθ̇2 = − µ
2
θ̈ + 2ṙθ̇ = 0
These are o
ital equations. For simplicity, use µ = 1.
Using same techniques as above come up with a system of equation in the form of
Ẏ = f(Y ). Here order your equations so Y = [r; ṙ; θ; θ̇].
Now start by writing an anonymous function O
itFunc(Y,mu) that takes as input the
vector Y and outputs f(Y ). Once you have this function we can now solve them using
two different techniques:
4
E7 Fall 2019, UC Berkeley
(a) Euler’s Method: Write a function Yout = Euler(Yinitial, t,∆t, O
itFunc) where
Yinitial is your initial conditions. Use Yinitial = [2; 1.5; 1.5; 3], ∆t = 0.001 and t=5.
Here Yinitial is the initial value of at t=0, t is the final time till which you need
to advance your solution to, ∆t is the time step size and O
itFunc returns f(Y).
Inside your function, you will step forward in time using Euler’s method. The
Euler’s equations, for value of Y at time tn+1, (where f(Y) is not a function of
time), is defined as:
Yn+1 = Yn + ∆tf(Yn)
where Yn is the Y value at the previous time step.
Define a vector tspan = 0 : ∆t : t. Your output Yout should be an a
ay of Y where
each column of Y is Y at that time step, i.e. number of columns of Y should be
equal to the length of tspan.
(b) RK4: Now lets use a much more stable time stepping scheme than Euler’s called
Runge-Kutta 4 (RK4) method. Again a write a function Yout = RK4(Yinitial, t,∆t,
O
itFunc). Use the same inputs as before and your Yout would be similar as well
as in the case of Euler’s. The method for time stepping will be changed however.
RK4 method for value of Y at time tn+1, (where f(Y) is not a function of time), is
defined as:
Yn+1 = Yn +
1
6
(k1 + 2k2 + 2k3 + k4)
where
k1 = ∆tf(Yn)
k2 = ∆tf(Yn + k1/2)
k3 = ∆tf(Yn + k2/2)
k4 = ∆tf(Yn + k3)
Store your calculated variables Y from both methods as YEulerO
and YRK4O
.
4. In this problem, you will use ode45, Euler and RK4 to simulate the o
its of a planet
around a (fixed) star and a moon around the planet, as illustrated in the figure below.
The equations of motion for this problem are quite complex, so we will use instead a
simplified set of equations to obtain a reasonable approximation of the system. The
dependent variables of these approximate equations of motion are: the star-planet dis-
tance R; the planet-moon distance r; the star-planet angle θ; and, the planet-moon
angle β. The simplified equations of motion are second-order and can be expressed as:
• R̈ = Rθ̇2 −GmSR−2 m·s−2 (radial acceleration of the planet)
5
E7 Fall 2019, UC Berkeley
Figure 2: Moon-Planet-Star O
it
• θ̈ = −2Ṙθ̇R−1 rad·s−2 (angular acceleration of the planet)
• r̈ = rβ̇2−GmP r−2−GmS(R+ r cos β)−2 m·s−2 (Radial acceleration of the moon)
• β̈ = −2ṙβ̇r−1 rad·s−2 (angular acceleration of the moon)
The constants that enter the preceding equations are:
• Mass of the star, mS = 1 × 1030 kg
• Mass of the planet, mP = 1 × 1023 kg
• Mass of the moon, mM = 1 × 1022 kg
• Gravitational constant, G = 6.673 × 10−11 m3·kg−1·s−2
Also, the initial conditions selected for the four variables are:
• Radial position of the planet about the star, rP/S = 1 × 1011 m
• Radial velocity of the planet, 0 m·s−1
• Angular position of the planet, 0 rad
• Angular velocity of the planet, ωP/S = 2.25 × 10−7 rad·s−1
• Radial position of the moon about the planet, rM/P = 1 × 109 m
• Radial velocity of the moon, 0 m·s−1
• Angular position of the moon, 0 rad
• Angular velocity of the moon, ωM/P = 2 × 10−6