AMATH 301 - Spring 2020
Homework #9
Due on Sunday, June 7, 2020
Instructions for submitting:
• Problems 1-5, 8 should be submitted to MATLAB Grader. You have 3 attempts
for each problem.
• Problems 6-7, 9-10 should be submitted to Gradescope. The solutions and the code
used to get those solutions should be submitted as a single pdf. All code should be
at the end of the file. You must select which page each problem is on when you
submit to Gradescope.
Linear Pendulum
Consider the linear second-order differential equation for the motion of a pendulum,
θ¨ = − g
L
θ.
Here θ (theta) is the angle of deflection of the pendulum from the vertical. g = 9.8
is acceleration due to gravity, L = 1 is the length of the pendulum. This differential
equation is a good approximation when θ is small.
By defining ω (omega) as the derivative ω = θ˙, we can convert this single second-order
ODE into a system of two first-order ODEs,
θ˙ = ω
ω˙ = − g
L
θ.
This system of ODEs is linear and therefore can be written in matrix form, θ˙
ω˙
 =
 0 1
− g
L
0
θ
ω
 .
Alternatively, we can write it as
x˙ = Ax
where
x =
θ
ω
 , A =
 0 1
− g
L
0
 .
We will solve the system of differential equations over the time interval 0 ≤ t ≤ T for
T = 10 with initial conditions θ(0) = 0, ω(0) = 0.5.
(5 points) Problem 1: MATLAB Grader
The forward Euler formula for a system of linear differential equations is
xk+1 = xk + ∆tAxk.
Implement the forward Euler method with ∆t = 0.005 to solve the initial value problem
for the pendulum. Make a 2 × 1 column vector with the forward Euler solution at the
final time T = 10.
The exact solution at time T = 10 is
x(10) =
θ(10)
ω(10)
 =
0.5√Lg sin (√ gL10)
0.5 cos
(√
g
L
10
)
 .
Calculate the exact solution. If the forward Euler solution at the final time is the vector
x FE and the exact solution at the final time is the vector x exact, then you can calculate
the global error of the forward Euler solution by doing
norm(x exact-x FE).
Calculate the global error and store it in the variable ans1.
(10 points) Problem 2: MATLAB Grader
The backward Euler formula for a system of linear differential equations is
xk+1 = xk + ∆tAxk+1.
Using some matrix algebra, we can write this as the following linear system
(I−∆tA)xk+1 = xk.
where I is the 2× 2 identity matrix. At each time step of backward Euler, we must solve
this linear system with the same matrix I −∆tA but a different right-hand side vector
xk. This is a good candidate for LU decomposition.
(a) Perform an LU decomposition on the matrix I − ∆tA with ∆t = 0.005 to obtain
matrices L, U, and P. Store the matrix L in the variable ans2.
(b) Implement the backward Euler method with ∆t = 0.005 to solve the initial value
problem for the pendulum. At each time step, you should use the LU decomposition
to solve the linear system. Calculate the global error and store it in the variable
ans3.
(5 points) Problem 3: MATLAB Grader
(a) Forward Euler will be stable for the pendulum problem if |λ| < 1 for all eigenvalues
of I+ ∆tA. Find the eigenvalue of I+ ∆tA that is largest in magnitude (or absolute
value) when ∆t = 0.005.
(b) Backward Euler will be stable for the pendulum problem if |λ| < 1 for all eigenvalues
of (I−∆tA)−1. Find the eigenvalue of (I−∆tA)−1 that is largest in magnitude (or
absolute value) when ∆t = 0.005. You may use the inv command to compute the
matrix inverse.
(c) Create a 1×2 row vector named ans4 that has the magnitude of the eigenvalue from
part (a) as the first component and the magnitude of the eigenvalue from part (b) as
the second component.
(10 points) Problem 4: MATLAB Grader
Recall that forward Euler is obtained by using a forward difference to approximate the
derivative and backward Euler is obtained by using a backward difference to approximate
the derivative. The Leapfrog method is based on using a central difference to approximate
the derivative:
xk+1 − xk−1
2∆t
= f(xk) ⇒ xk+1 = xk−1 + 2∆tf(xk).
This method is explicit like forward Euler because the formula we end up with only
depends on past solution values. A major difference between this method and the Euler
methods is that it is a multistep method – the iteration actually uses two past solution
values: xk−1 and xk. This means that each time we generate a new iterate xk+1, we need
to use the ‘current’ solution value xk as well as the one before it, xk−1.
Since this method is a multistep method, we cannot start the iteration with just the
initial values because the formula to compute x1 is
x1 = x−1 + 2∆tf(x0),
which is a problem since x−1 isn’t defined. To get around this, we perform a single step
of Forward Euler to calculate x1, and then switch over to Leapfrog to calculate x2 and
the remaining iterates.
Implement the Leapfrog method with ∆t = 0.005 to solve the initial value problem for
the pendulum. Make a 2 × 1 column vector named ans5 with the Leapfrog solution at
the final time T = 10. Calculate the global error and store it in the variable ans6.
(10 points) Problem 5: MATLAB Grader
Use the MATLAB function ode45 to solve the initial value problem for the pendulum.
Make a 1× 2 row vector named ans7 with the ode45 solution at the final time T = 10.
Calculate the global error and store it in the variable ans8.
Hint: When you subtract two vectors to calculate the error, you must have either two
row vectors or two column vectors. You cannot have one of each.
(10 points) Problem 6: Gradescope
Construct a phase portrait with θ on the horizontal axis and ω on the vertical axis. Use
meshgrid to generate a grid of points between −0.75 ≤ θ ≤ 0.75 and −0.75 ≤ ω ≤ 0.75
with 25 equally-spaced points in both directions. Use the quiver function to draw a grid
of arrows with components (θ˙, ω˙), where θ˙ and ω˙ are given by the ODEs of the pendulum
system. Include labels for the axes. Unlike the activity, you do not need to draw axes.
You should see something (somewhat) like the following:
-0.6 -0.4 -0.2 0 0.2 0.4 0.6
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
The exact solution to the system of differential equations with initial conditions
θ(0) = 0 and ω(0) = 0.5 is
θ(t) = 0.5
√
L
g
sin
(√
g
L
t
)
ω(t) = 0.5 cos
(√
g
L
t
)
Create a vector of the exact solution θ(t) for all times from t = 0 to t = 10 in steps of
∆t = 0.005. Do the same for ω(t). Use these vectors to add a trajectory to the phase
portrait. Make the color of the trajectory different from the color of the arrows.
(15 points) Problem 7: Gradescope
We will now construct another phase portrait, but instead of drawing a trajectory for the
exact solution we will draw trajectories using numerical solutions.
(a) Begin exactly as in Problem 6. Use meshgrid and quiver to draw the same grid of
arrows. Then label the axes.
(b) Implement the forward Euler method with ∆t = 0.005 to solve the initial value prob-
lem for the pendulum. Create vectors that contain the forward Euler approximations
of θ(t) and ω(t) at all time steps from t = 0 to t = 10. Use these vectors to add a
trajectory to the phase portrait.
(c) Repeat part (b) for both the backward Euler and Leapfrog methods.
(d) All trajectories should be different colors. Include a legend so that the methods can
be identified. In order to include the trajectories in your legend and not the quiver
arrows, you may want to look at the documentation for the legend command in the
section “Included Subset of Graphics Objects in Legend”.
Nonlinear Pendulum
The linear pendulum equation is based on the small angle approximation sin(θ) ≈ θ.
Without this approximation, the equation is more accurate and becomes
θ¨ = − g
L
sin(θ).
This second-order ODE is nonlinear and so is more difficult to solve. In particular,
implicit methods like backward Euler are harder to implement for nonlinear ODEs.
This equation can also be converted to a system of differential equations in terms of θ and
ω = θ˙. We will solve the system of differential equations for the nonlinear pendulum over
the same time interval 0 ≤ t ≤ T for T = 10 with the same initial conditions θ(0) = 0,
ω(0) = 0.5.
(10 points) Problem 8: MATLAB Grader
Use the MATLAB function ode45 to solve the initial value problem for the nonlinear
pendulum. Make it so that ode45 outputs the solutions at equally spaced times from
t = 0 to t = 10 in steps of ∆t = 0.005. Create a 2001 × 1 column vector named ans9
that contains the values of θ(t) at each time. Then create a column vector named ans10
that contains the values of ω(t).
(20 points) Problem 9: Gradescope
Construct a phase portrait for the nonlinear pendulum. Use meshgrid to generate a grid
of points. Use 25 equally-spaced points for θ between −2pi and 2pi. Use 20 equally-spaced
points between −8 and 8 for ω. Use the quiver function to draw a grid of arrows. The
resulting picture should look similar to what you obtained in Problem 6 in the middle of
the plot, but different around the edges (this is due to the nonlinearity of the system).
Use ode45 to solve the system at the same times as in Problem 8 and with the following
variety of initial conditions:
• θ(0) = 0, ω(0) = 0.5
• θ(0) = 2, ω(0) = 1
• θ(0) = pi, ω(0) = −10−4
• θ(0) = 2pi, ω(0) = −7
• θ(0) = −2pi, ω(0) = 7
Add these solution trajectories to the phase portrait. These solution trajectories should
flow with the arrows from quiver. Set the axes to display from −2pi < θ < 2pi and
−8 < ω < 8. You do not need to add a legend because it is self-explanatory which
trajectory belongs to which initial condition.
These trajectories come in two varieties: those that repeatedly encircle the origin, and
those that do not. Describe the behavior of the pendulum for each of these two types of
trajectories. Remember that θ represents an angle so θ = 0 is the same as θ = 2pi and
θ = pi is the same as θ = −pi.
(5 points) Problem 10: Gradescope
Which of the following methods are explicit methods? Choose all that apply.
(a) forward Euler
(b) backward Euler
(c) midpoint method (RK2)
(d) fourth-order Runge-Kutta (RK4)
Gradescope Deliverables Your Gradescope writeup should contain the following:
• Problem 6: The phase portrait
• Problem 7: The phase portrait
• Problem 9: The phase portrait and a description of the behavior for each type of
trajectory
• Problem 10: A letter or multiple letters corresponding to your answer choice
• Code: Code for problems 6, 7 and 9