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作业君