-------------------------------------------------------------------------
PHYS210 2009-11-05 LAB SUMMARY
1) Discussion of the critical value, omega0*, for the nonlinear
pendulum
2) Modification of pendulum.m (-> pendulume.m) to compute energies
3) Convergence test of total energy conservation of nonlinear pendulum
4) Free time to work on homework, term projects
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1) Discussion of the critical value, omega0, for the nonlinear pendulum
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Here are the critical values omega0* that I computed for various levels
-------------------------------------------------------------------------
Critical omega0 value as function of level
-------------------------------------------------------------------------
level omega0* d(omega0*) 4^(level - 6) d(omega0*)
6 1.929879441857338 0.053608 0.0536
7 1.983487772196532 0.0124290 0.0497
8 1.995916748046875 0.0030649 0.0490
9 1.998981679230927 7.63802e-04 0.0489
10 1.999745482113212 1.90935e-04 0.0489
11 1.999936417131920
-------------------------------------------------------------------------
where d(omega0*) = omega0*(level+1) - omega0*(level)
There are two things that are apparent from this table
1) omega0* is converging to some value at a rate which is O(delt^2)
2) It is natural to conjecture that the continuum value of omega0*
is 2
Indeed, from today's lecture discussion of the various energy quantities
associated with the pendulum we have
T = omega^2/2
V = 1 - cos(theta)
E = T + V
for the kinetic, potential and total energy, respectively, and where
we recall that we are working in units in which m = g = L = 1.
Thus for the critical case, where we give the pendulum the precise initial
kick needed for it to swing around by an angle pi and stop there (in
principle, forever, even though it is an unstable equilibrium), we have
E_initial = T_initial + V_initial
= (omega0*)^2/2 + 0 = (omega0*)^2/2
E_final = T_final + V_final
= 0 + (1 - cos(pi)) = 2
Now, since E_initial = E_final we have
(omega0*)^2/2 = 2 -> omega0* = 2
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2) Modification of pendulum.m (-> pendulume.m) to compute energies
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PART 1: CREATING pendulume.m
Copy pendulum.m from ~phys210/pendulum to ~/octave/pendulume.m
NOTE the change in the name of the function file---pendulume.m (with an
'e' for energy after pendulum), rather than pendulum.e
% cd
% cp ~phys210/octave/pendulum.m ~/octave/pendulume.m
Modify pendulume.m so that it defines the function pendulume with the
following header
function [t theta omega T V E] = pendulume(tmax, level, theta0, omega0)
where the additional output arguments, relative to pendulum, are T, V and E,
and are defined as follows
% T: (real vector) Vector of length nt containing computed
% kinetic energy at discrete times t(n).
% V: (real vector) Vector of length nt containing computed
% potential energy at discrete times t(n).
% E: (real vector) Vector of length nt containing computed
% total energy at discrete times t(n).
Add code to the function definition to compute T, V and E using the
formulae derived in class (for m = g = L = 1) and already mentioned
above
T = omega^2/2
V = 1 - cos(theta)
E = T + V
Hint: You can compute all of the above using whole-array calculations
(using element-by-element operators as necessary) after the for loop and
the definition of omega(nt). For completeness, you should also change
the occurrences of 'pendulum' in the fprintf statement to 'pendulume'
PART 2: BASIC TEST OF pendulume
Create a script file ~/octave/tpe.m (i.e. ensure that tpe.m is created
in your ~/octave directory) that runs pendulume with the following
parameters
tmax = 40.0
theta0 = 0
omega0 = 1.9
level = 6
and then graphs T, V and E vs t on the same plot with T in red,
V in green and E in blue
Recall some basic plotting commands
clf; Clears the current figure window
hold on; Do NO clear the figure window each time a new plot is made
plot(x,y,'r') Plot the vector y vs the vector x using a red line
Available colors are
'r' red
'g' green
'b' blue
'c' cyan
'm' magenta
'k' black
'w' white
Making a script pause and wait for user input
input('Type ENTER to continue: ');
Run the script by typing
>> tpe
at the prompt and modify tpe.m and pendulume.m as necessary until
your plot looks reasonable.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3) Convergence test of total energy conservation of nonlinear pendulum
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PART 1: Investigating convergence behaviour of deviation of total
energy
Add code to tpe.m so that it runs pendulume with the same basic input
parameters
tmax = 40.0
theta0 = 0
omega0 = 1.9
but for levels 6, 7, 8 and 9
After each invocation of pendulume, compute a vector containing the
deviation in the total energy, relative to the initial total energy,
as a function of the discrete time using a statement such as
dE = E - E(1)
Note that dE will then be a vector of length nt
Have tpe.m make a single plot showing all of these deviations, using a
different colored line for each dataset.
NOTE: GIVEN THAT YOU NEED TO MAKE A SINGLE PLOT WITH DIFFERENT COLORED
LINES, IT IS MOST STRAIGHTFORWARD TO DO THE CONVERGENCE TEST *WITHOUT*
USING A for loop.
What do you observe?
PART 2: Investigating convergence behaviour of deviation of total
energy by scaling deviations
Add additional code to tpe.m so that it performs the same test
described above but so that it now computes a SCALED deviation of
the total energy using
dEs = 4^(level - 6) * (E - E(1))
for level = 6, 7, 8 and 9
Again, make a single plot showing all of these deviations, using
different colors for different datasets.
What do you observe?