Programming Assignment 2: Numerical Convergence

ESM 3124 - Fall 2020

Due: 11-13-2020

Two results of your first programming assignment are functions that can be called with:

dSdt = my_ode(t, S, g, L)

[S, t] = eulers_method(@(t, S) my_ode(t, S, g, L),tspan, S0, delta_t)

Several students observed that the total energy of the system, calculated using their numerical

solution, grew with time. This apparent violation of the conservation of energy has nothing to do

with the pendulum being studied, and instead depends on the details of the numerical solution.

Notice that eulers_method(...,delta_t) is a function of delta_t, the time step chosen to

simulate the motion. This assignment will focus on determining exactly how the solution depends

on delta_t. The assignment is broken into 2 main parts. All submissions need to be in a single

.zip folder labeled 'Your_Name_Pendulum_Convergence.zip'.

0.1 Convergence Analysis (10 pts)

A necessary condition for calculating an accurate numerical solution is that your results are

convergent. Conceptually, this means that as the time step ∆t used to simulate motion becomes

smaller, we expect that the solution converges to a single state. If we have two solutions θi(t)

and θi+1(t) that are calculated using time steps ∆ti and ∆ti+1, we can calculate the magnitude

of the difference in the solutions at a fixed point in time ∆θi = |θi+1(t0)− θi(t0)|. If the solutions

converge to a single value, then ∆θi → 0 as ∆ti → 0. A solution satisfying this property is said

to be convergent. Now that you've finished reading this paragraph re-read this paragraph, it

explains the entire assignment.

While this may make us want to choose the smallest possible ∆t, accuracy comes with the trade-off

of increased run time. Your goal in this section is to design and write a function that determines

if your solution is convergent, and times how long it takes eulers_method() to run for numerous

values of ∆t. Assess convergence using 2 fairly simple metrics 1) the final angle θf (t = tf ), and

the final system energy Ef (t = tf ). Use your code from the first programming assignment as

the starting point for this, but remove any parts that are not used in this assignment (e.g., the

animation section).

1. Write a local function:

[del_S, del_E, run_time] = convergence_analysis(..., delta_t_array)

that measures the convergence of your eulers_method function by looking at the change

in the final angle and energy as a function of delta_t. Your function must also measure

run_time, how long it takes to simulate each trajectory. (See the help menu for tic and toc

to learn how to time your code)

* You will need to decide what most of the function inputs will be. However, one input argu-

ment to convergence_analysis() must be a columnn vector delta_t_array that contains

all values of ∆ti that you will simulate motion for.

0.2 Determine Valid ∆t (10 pts)

Use the function written above to simulate trajectories for 30 logarithmically spaced values from

∆t = 10−5 seconds to ∆t = 10−1 seconds (pull up the help menu for logspace() for assistance).

Plot how the solution converges as a function of ∆t, and what the associated run_times are. Using

your plots, conclude what values of ∆t are useful for simulating motion.

1. For each value of ∆t, simulate 10 seconds of pendulum motion, then save the final values

of the angle theta_f and total energy E_f. After calculating these values for each delta_t

in delta_t_array, calculate the change in each value as the delta_t increases. Simulate a

pendulum with L = 1m, M = 1kg, and initial conditions θ0 = pi/1.1 and θ˙0 = 0.

1

Figure 1: Example of convergence results

2. Plot the result of your calculations in a single figure with two subplots. One plot must show

the convergence of the final angle ∆θf vs ∆t, and the second the final energy ∆Ef vs ∆t.

(See loglog() for plotting logarithmically spaced values, and yyaxis for creating different

left / right vertical axes).

3. Print a 200 dpi .png image of your plot called convergence_plot.png.

Submission Format

Make sure your code is properly commented to explain what each portion does, and how.

In a single .zip folder named `Your_Name_Pendulum_Convergence.zip' you will submit:

1. A single .m file that contains your code and local functions described in Part 0.1. The TA

and Professor Domann should be able to load this code into Matlab, press Run, and have

your code execute successfully.

2. An image called convergence_plot.png that contains an image of your convergence analysis.

2

Recommended Code Update:

Compare how you used ode23() and eulers_method() in the first programming assignment.

odefun = @(t,S) my_ode(t,S,g,L);

[t, S] = ode23(odefun, tspan, S0, options);

[t, S] = eulers_method(odefun, tspan, S0, delta_t);

These functions are structured to do the same thing, numerically simulate an ODE. Note that

the inputs to both functions are nearly identical. The only information ode23() receives about

the ODE it is simulating is odefun = @(t, S) my_ode(t, S, g, L), and the initial state S0. In

other words, ode23() is flexible enough to make calculations using only those properties. Now

consider how most people wrote eulers_method()

function [t, S] = eulers_method(odefun, tspan, S0, delta_t)|

%% Function Setup

...

%% Run Problem Specific Loop

for i = 1 : numel(t)-1

% rate of change of state

dsdt = odefun(t(i), S(i,:));

% update angle and angular velocity

theta(i+1) = theta(i) + dsdt(1)*delta_t;

omega(i+1) = omega(i) + dsdt(2)*delta_t;

% update state

S(i+1,1) = theta(i+1);

S(i+1,2) = omega(i+1);

end %end loop

end %end function

While that code is very readable, it only works for a problem with 1 position variable (θ) and one

velocity variable (ω). What if you want to simulate a 2 degree of freedom system like a double

pendulum? Your code would no longer work. Fortunately there is a very easy change. Note that

the code for % update angle ... merely serves as temporary storage, and is not actually needed.

As a result the code can simultaneously be condensed, and made more flexible.

function [t, S] = eulers_method(odefun, tspan, S0, delta_t)|

%% Function Setup

...

%% Run More Flexible Loop

for i = 1 : numel(t)-1

% rate of change of state

dsdt = odefun(t(i), S(i,:));

% update state

S(i+1,:) = S(i, :) + dsdt * delta_t;

end %end loop

end %end function

As long as you initialize S so it has the same number of columns that dsdt has elements, then

this code will work. Note that the code S(i, :) is read `S at row i, all columns'. If you make

this change, then when you want to simulate a new system the only code you need to change is

the definition of odefun(), where you write your system of 1st order ODEs. In other words, your

numerical method eulers_method() no longer depends on a single problem, and is instead a tool

capable of simulating any properly structured ODE. I recommend doing this now before the third

and final programming assignment.

3

欢迎咨询51作业君