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