pendulum and lpendulumStart a Matlab session and arrange the command window so that you can work in it while viewing these notes, since you'll be doing a lot of cut and pasting in the early part of the lab. We won't need a terminal window for the most part.
BROWSE SOURCE
pendulum.m (Nonlinear)lpendulum.m (Linear)pendulum solves the nonlinear pendulum equation using
\(O(h^2)\) finite difference techniques as discussed in class.
As usual, we can examine the function
definition from the Matlab command line using type:
>> type pendulum
function [t theta omega] = pendulum(tmax, level, theta0, omega0)
% Solves nonlinear pendulum equation using second order finite difference
% approximation as discussed in class.
%
% Input arguments
%
% tmax: (real scalar) Final solution time.
% level: (integer scalar) Discretization level.
% theta0: (real scalar) Initial angular displacement of pendulum.
% omega0: (real scalar) Initial angular velocity of pendulum.
%
% Return values
%
% t: (real vector) Vector of length nt = 2^level + 1 containing
% discrete times (time mesh).
% theta: (real vector) Vector of length nt containing computed
% angular displacement at discrete times t(n).
% omega: (real vector) Vector of length nt containing computed
% angular velocity at discrete times t(n).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% trace controls "tracing" output. Set 0 to disable, non-0 to enable.
trace = 1;
% tracefreq controls frequency of tracing output in main time step loop.
tracefreq = 100;
if trace
fprintf('In pendulum: Argument dump follows\n');
tmax, level, theta0, omega0
end
% Define number of time steps and create t, theta and omega arrays of
% appropriate size for efficiency.
nt = 2^level + 1;
t = linspace(0.0, tmax, nt);
theta = zeros(1,nt);
omega = zeros(1,nt);
% Determine discrete time step from t array.
deltat = t(2) - t(1);
% Initialize first two values of the pendulum's angular displacement.
theta(1) = theta0;
theta(2) = theta0 + deltat * omega0 - 0.5 * deltat^2 * sin(theta0);
if trace
fprintf('deltat=%g theta(1)=%g theta(2)=%g\n',...
deltat, theta(1), theta(2));
end
% Initialize first value of the angular velocity.
omega(1) = omega0;
% Evolve the oscillator to the final time using the discrete equations
% of motion. Also compute an estimate of the angular velocity at
% each time step.
for n = 2 : nt-1
% This generates tracing output every 'tracefreq' steps.
if rem(n, tracefreq) == 0
fprintf('pendulum: Step %d of %d\n', n, nt);
end
theta(n+1) = 2 * theta(n) - theta(n-1) - deltat^2 * sin(theta(n));
omega(n) = (theta(n+1) - theta(n-1)) / (2 * deltat);
end
% Use linear extrapolation to determine the value of omega at the
% final time step.
omega(nt) = 2 * omega(nt-1) - omega(nt-2);
end
lpendulum is a precisely analogous function that solves the
linear pendulum equation (the one you've seen in previous
physics courses) and, apart from comments, diagnostic
messages, and variable name changes, the code for
the two functions differ in literally two lines.
>> type lpendulum
function [t theta omega] = lpendulum(tmax, level, theta0, omega0)
% Solves linear pendulum equation using second order finite difference
% approximation.
%
% Input arguments
%
% tmax: (real scalar) Final solution time.
% level: (integer scalar) Discretization level.
% theta0: (real scalar) Initial angular displacement of pendulum.
% omega0: (real scalar) Initial angular velocity of pendulum.
%
% Return values
%
% t: (real vector) Vector of length nt = 2^level + 1 containing
% discrete times (time mesh).
% theta: (real vector) Vector of length nt containing computed
% angular displacement at discrete times t(n).
% omega: (real vector) Vector of length nt containing computed
% angular velocity at discrete times t(n).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% trace controls "tracing" output. Set 0 to disable, non-0 to enable.
trace = 1;
% tracefreq controls frequency of tracing output in main time step loop.
tracefreq = 100;
if trace
fprintf('In lpendulum: Argument dump follows\n');
tmax, level, theta0, omega0
end
% Define number of time steps and create t, theta and omega arrays of
% appropriate size for efficiency.
nt = 2^level + 1;
t = linspace(0.0, tmax, nt);
theta = zeros(1, nt);
omega = zeros(1, nt);
% Determine discrete time step from t array.
deltat = t(2) - t(1);
% Initialize first two values of the pendulum's angular displacement.
theta(1) = theta0;
theta(2) = theta0 + deltat * omega0 - 0.5 * deltat^2 * theta0;
if trace
fprintf('deltat=%g theta(1)=%g theta(2)=%g\n',...
deltat, theta(1), theta(2));
end
% Initialize first value of the angular velocity.
omega(1) = omega0;
% Evolve the oscillator to the final time using the discrete equations
% of motion. Also compute an estimate of the angular velocity at
% each time step.
for n = 2 : nt-1
% This generates tracing output every 'tracefreq' steps.
if rem(n, tracefreq) == 0
fprintf('lpendulum: Step %d of %d\n', n, nt);
end
% Note that the last term is deltat^2 * theta(n), whereas
% for the nonlinear pendulum it is deltat^2 * sin(theta(n)).
% This is the ONLY difference in the FDAs for the two
% systems.
theta(n+1) = 2 * theta(n) - theta(n-1) - deltat^2 * theta(n);
omega(n) = (theta(n+1) - theta(n-1)) / (2 * deltat);
end
% Use linear extrapolation to determine the value of omega at the
% final time step.
omega(nt) = 2 * omega(nt-1) - omega(nt-2);
end
We will first use these functions in tandem so that we can compare the behaviour of the pendulum for the nonlinear/linear cases.
In all of the calculations below, we will start with
theta0 = 0 (i.e. so that the pendulum is hanging vertically),
but will vary the initial angular velocity omega0
Assign values to variables that will be used as input arguments to the functions---be sure to type all of the assignment statements carefully!
Start with a small value of omega0
>> tmax = 40.0
>> level = 8
>> theta0 = 0
>> omega0 = 0.01
Recall that level implicitly defines
via
\[ n_t = 2^{\rm level + 1} \]
\[ \Delta t = \frac{t_{\rm max}}{2^{\rm level}} \]
so in this case, the solution will be computed at a total of 257 discrete times, including the first and second time steps where the values are determined by the initial conditions.
Now call pendulum, being careful to specify a vector of three names on the left hand side of the assignment so that all output arguments from pendulum are captured
>> [t theta omega] = pendulum(tmax, level, theta0, omega0);
In pendulum: Argument dump follows
tmax = 40
level = 8
theta0 = 0
omega0 = 0.010000
deltat=0.15625 theta(1)=0 theta(2)=0.0015625
pendulum: Step 100 of 257
pendulum: Step 200 of 257
Now call lpendulum with the same arguments, but being sure to
use different names on the left hand side for the angular
displacement and velocity , i.e. ltheta and lomega. It
doesn't matter that we use the same name t in the two calls
since the output vector of times is the same for the two calculations.
>> [t ltheta lomega] = lpendulum(tmax, level, theta0, omega0);
In lpendulum: Argument dump follows
tmax = 40
level = 8
theta0 = 0
omega0 = 0.010000
deltat=0.15625 theta(1)=0 theta(2)=0.0015625
lpendulum: Step 100 of 257
lpendulum: Step 200 of 257
Plot the results: note the use of two plot commands with
hold on: recall that the latter tells Matlab not to clear the figure
window each time a plot command is executed
>> clf; hold on; plot(t, theta, 'r'); plot(t, ltheta, 'og');
Observe that the period of the oscillation is \(2\pi\), which follows from our use of "natural units", \(g=L=1\).
Also note how there is very little difference in the motions of the two types of pendula; this is expected since the amplitude is small in this case---about 0.01 radians---so \(\theta\) and \(\sin(\theta)\) are nearly equal at all times
With Matlab's ability to perform whole-array operations, it is easy to calculate and plot the actual difference in the \(\theta\)'s between the two simulations.
>> clf; plot(t, theta-ltheta);
Note that the difference is oscillatory, with an amplitude that increases in time---the difference is real, i.e. it is not due to the approximate nature of the solutions but represents the small deviation of the nonlinear pendulum from precisely-linear behaviour.
CONVERGENCE TESTING!
As a stringent test that our numerical solution is behaving
as it should---implying that the implementation is correct---we
will now perform a basic 3-level convergence test for the nonlinear
pendulum, as discussed in class. That is, we will perform a set of
three calculations where three of the four parameters that we supply
to pendulum, namely
tmaxtheta0omega0are held fixed, while the grid spacing (time step),
\(\Delta t\),
is varied. Again,
keep in mind that we change the time step (by factors of two)
by changing the integer-valued level parameter.
For this test we will use
From this point onward, I recommend that you use "copy and paste" to enter the commands, and to expedite this I will include he prompt at the beginning of each Matlab command line, This will allow you to select and copy an entire block of commands from your browser, and then paste them into the Matlab session. Just be careful not to copy any of the Matlab output from the command/output block.
However, if you do choose to type the commands manually, be careful with your typing, paying particular attention to the names on the left hand sides of the assignment statements
[t6 theta6 omega6] = pendulum(tmax, 6, theta0, omega0);
In pendulum: Argument dump follows
tmax = 40
level = 6
theta0 = 0
omega0 = 0.010000
deltat=0.625 theta(1)=0 theta(2)=0.00625
[t7 theta7 omega7] = pendulum(tmax, 7, theta0, omega0);
In pendulum: Argument dump follows
tmax = 40
level = 7
theta0 = 0
omega0 = 0.010000
deltat=0.3125 theta(1)=0 theta(2)=0.003125
pendulum: Step 100 of 129
[t8 theta8 omega8] = pendulum(tmax, 8, theta0, omega0);
In pendulum: Argument dump follows
tmax = 40
level = 8
theta0 = 0
omega0 = 0.010000
deltat=0.15625 theta(1)=0 theta(2)=0.0015625
pendulum: Step 100 of 257
pendulum: Step 200 of 257
First, graph \(\theta(t)\) from the 3 calculations on the same plot:
clf; hold on; plot(t6, theta6, 'r-.o');
plot(t7, theta7, 'g-.+'); plot(t8, theta8, 'b-.*');
Note that the curves show general agreement, but that as time progresses there is an increasing amount of deviation, which is most pronounced for the calculation that used the largest grid spacing (level 6, red line)
The observed deviations are a result of the approximate nature of the finite difference scheme
Now, "downsample" the level 7 and 8 theta arrays so that they have the
same length as the level 6 arrays (with values defined at the
times in t6). This means we select every second element of theta7
and every fourth of theta8, which is easy to do using the colon
operator:
theta7 = theta7(1:2:end);
theta8 = theta8(1:4:end);
Now, compute the differences in the grid functions from level to level, assigning them to new variables ...
dtheta67 = theta6 - theta7;
dtheta78 = theta7 - theta8;
and make a plot of the differences.
clf; hold on; plot(t6, dtheta67, 'r-.o'); plot(t6, dtheta78, 'g-.+');
Note that the two curves have the same shape, but different overall amplitudes; each curve is (up to some numeric factor that is approximately 1, i.e. "of order unity") a direct measure of the error in the numerical solution
Now scale (multiply) dtheta78 by 4 and replot together
with dtheta67. If the convergence is second order, the curves should be
nearly coincident (the grid spacing goes down by a factor of 2,
so the differences should go down by a factor of \(2^2 = 4\))
dtheta78 = 4 * dtheta78;
clf; hold on; plot(t6, dtheta67, 'r-.o'); plot(t6, dtheta78, 'g-.+');
Are they nearly coincident?
YES!!
Also note how the magnitude of the error increases with time, which is characteristic of a finite difference solution of this type
We can extend our convergence test by repeating the calculation
with an even higher level, level = 9.
Again, recall that this means the grid spacing---the discrete time step, \(\Delta t\)---will be reduced by a factor of 2 relative to the level 8 computation.
level = 9
[t9 theta9 omega9] = pendulum(tmax, 9, theta0, omega0);
Downsample theta9, compute the difference relative to theta8, scale by the appropriate power of 4, which is \(4^2 = 16\), and add to the plot:
theta9 = theta9(1:8:length(theta9));
dtheta89 = theta8 - theta9;
dtheta89 = 16 * dtheta89;
plot(t6, dtheta89, 'b-.*');
Note how the sequence from red -> green -> blue curves shows
increasingly good overlap: this is exactly what we expect, since
as \(\Delta t \to 0\),
i.e. as level \(\to \infty\), the leading
order error term in the finite difference solution
which is \(O(\Delta t^2)\),
will become increasingly dominant over the higher order terms
\(O(\Delta t^4), O(\Delta t^6). \ldots\)
We could continue the convergence test to even higher levels,
level = 10, 11, ... , but we should note that for sufficiently
high level (sufficiently small \(\Delta t\)) the finite difference
approximation would not improve, and, in fact, would begin to
degrade as a result of "round off" effects stemming from the
fact that the numbers being manipulated have a finite precision
(finite number of digits). This is an advanced topic that we
will not dwell on here.
In any case, from the convergence test, we now have strong evidence that the implementation is correct
We can therefore move on to "doing some science" with the code; ultimately, of course, this is the purpose of any project in computational science
In other words, we now will carry out some "numerical experiments" to investigate/compare the dynamical behaviour of the pendula.
We start by performing simulations using a sequence of increasing
values of omega0 , for fixed tmax and
theta0, and at level = 8.
Physically, this means that we will be giving the pendula ever increasing initial angular velocities, so the resulting amplitude of the motion will also increase
A key result that you can see immediately is whereas the period (frequency) of the linear pendulum (green) stays constant at \(2\pi\), i.e. is not dependent on the amplitude of the oscillation, the period of the nonlinear pendulum (red) does depend on the amplitude (or, equivalently, on omega0)
Also note that in addition to plots of theta, we will also be generating "phase space" graphs, in which \(\omega(t)\) is plotted against \(\theta(t)\) (with \(t\) being a parametric independent variable). Phase space analysis is extensively used in the study of dynamical systems.
Again, you are urged to use copy-and-paste to enter the commands, since we will be working through the notes fairly quickly:
First reset the fixed parameters, just as a precaution
tmax = 40; level = 8; theta0 = 0;
Now, use a loop to perform simulations of both the nonlinear and linear pendula, for various values of omega0 --- be sure to give yourself some time to study each plot before you press enter in response to the prompt
Copy and paste this entire block of statements:
for omega0 = [0.1 0.3 1.0 1.5 1.9]
[t theta omega] = pendulum(tmax, level, theta0, omega0);
[t ltheta lomega] = lpendulum(tmax, level, theta0, omega0);
clf; hold on; ptitle = sprintf('omega_0 = %g',omega0);
plot(t, theta, 'r'); plot(t, ltheta, '-.og');
title(ptitle);
input('Type ENTER to continue: ');
clf; hold on; ptitle = sprintf('Phase space plot: omega_0 = %g',omega0);
plot(theta, omega, 'r'); plot(ltheta, lomega, '-.og');
title(ptitle);
xlabel('theta');
ylabel('omega');
input('Type ENTER to continue: ');
end
Now, let's give the nonlinear pendulum an even larger initial kick ..
omega0 = 2.05;
[t theta omega] = pendulum(tmax, level, theta0, omega0);
clf; plot(t,theta,'r');
What is the physical interpretation of the behaviour of \(\theta(t)\) in this case?
Note that we have now seen two distinct types of behaviour in the nonlinear oscillator as illustrated by the following ...
[t thetalo omegalo] = pendulum(tmax, level, theta0, 1.9);
[t thetahi omegahi] = pendulum(tmax, level, theta0, 2.05);
clf; hold on;
plot(t, thetalo, 'g'); plot(t, thetahi, 'r');
title('"low" and "high" behaviour');
so somewhere in the interval [1.9, 2.05] there should be a "critical value" of omega0 where the behaviour type changes.
We will try to find this value using the technique of binary search, also known as bisection.
ASIDE
As mentioned above, most of the calculations that we
have performed thus far are encoded in a script tp that you
should be able to execute via
>> tp
You may want to examine the contents of tp.m using type tp
or browse it here
as an another example of a basic driver script (main program) that calls a function to do the core calculations of one or more simulation, and which then analyzes, plots, etc. the results.
END ASIDE
There are a set of course script and functions, defined
in ~/phys210/matlab that can be used to expedite a binary search
(google "bisection method" if you are still confused about this topic after
our discussion of it).
To see their names and synopses, use the lookfor command
>> lookfor phys210-bs
bsbracket - Print current binary search bracket [PHYS210-bs]
bscurrent - Return the midpoint of binary search bracket [PHYS210-bs]
bsdisp - Print values used in binary search [PHYS210-bs]
bsdone - Returns 1 [0] if binary search is done [not done] [PHYS210-bs]
bsfraction - Returns fractional width of binary search bracket [PHYS210-bs]
bshi - Replaces high value of binary search bracket with midpoint [PHYS210-bs]
bshigh - bshi Replaces high value of binary search bracket with midpoint [PHYS210-bs]
bsinitialized - Returns 1 [0] if binary search has [has not] been initialized [PHYS210-bs]
bslo - Replaces low value of binary search bracket with midpoint [PHYS210-bs]
bslow - bslo Replaces low value of binary search bracket with midpoint [PHYS210-bs]
bsnew - Initializes binary search [PHYS210-bs]
bsreset - Reset binary search [PHYS210-bs]
bssettol - Set tolerance for binary search [PHYS210-bs]
bstol - Return tolerance for binary search [PHYS210-bs]
bsundo - Restores previous midpoint of search [PHYS210-bs]
I've ended all of the H1 lines in the initial comment blocks in the
function definitions with [PHYS210-bs] so that the lookfor phys210-bs
will match all of them (the lookfor search is case-insensitive).
The ones that we will be using are
bsnewbslow (has an alias bslo)bshigh (has an alias bshi)bsundo (for undoing mistakes made when doing a searchbscurrentand less frequently
bsbracketbsfractionand they are invoked and operate as follows:
bsnew(<low>, <high>)
Initializes a new binary search in the current directory with an initial bracket [<lo>, <hi>]
bslow
Updates bracket so that the current midpoint becomes the new low value.
bshigh
Updates bracket so that the current midpoint becomes the new high value.
bsundo
Restores the previous bracket, so if you entered bslow when you
should have entered bshigh, type bsundo followed by `bshigh` and
vice versa. Can be applied repeatedly to return to any bracket
generated during the search.
bscurrent
Returns the current midpoint of bracket
bsbracket
Displays current endpoint values (low and high)
bsfraction
Returns normalized width of current bracket. Effectively tells you
to how many digits you have narrowed the search. When bsfraction
returns something of the order of 1.0e-16 (i.e. roughly `eps`),
you have narrowed the search about as far as you can
Type all of the following commands within a Matlab session.
Initialize the search: Note that the script returns the values defining the initial interval
bsnew(1.9, 2.05)
1.8999999999999999e+00
2.0499999999999998e+00
1.9749999999999999e+00
The function displays the two values that define the initial bracket, then an empty line and then the midpoint of the bracket
Let's first ensure that the fixed problem parameters are set: in
particular, we will perform the search using level = 7
tmax = 40; level = 7; theta0 = 0;
Set omega0 to the midpoint of the interval
omega0 = 1.9749999999999999e+00
Solve the pendulum equation, plot the angular displacement and, from the plot, determine which type of behaviour it exhibited
[t theta omega] = pendulum(tmax, level, theta0, omega0); clf; plot(t,theta);
In this case you should see the displacement oscillating---i.e.
"low" behaviour---so execute the function bslow to replace the
previous lower value for the bracket with the midpoint, and return
the new midpoint;
bslow
1.9749999999999999e+00
2.0499999999999998e+00
2.0124999999999997e+00
As with bsnew, bslow displays the current (updated) bracket and
the new midpoint value.
Use that new value in the next simulation:
omega0 = 2.0124999999999997e+00
[t theta omega] = pendulum(tmax, level, theta0, omega0); clf; plot(t,theta);
This time you should observe "high" behaviour---i.e. theta
monotonically increasing with time, so this time we call the function
bshigh:
bshigh
1.9749999999999999e+00
2.0124999999999997e+00
1.9937499999999999e+00
This function also displays the updated bracket (note that it is the high endpoint that has been changed this time), and the new midpoint.
Keep going ... set omega0 to the new midpoint and simulate:
omega0 = 1.9937499999999999e+00
[t theta omega] = pendulum(tmax, level, theta0, omega0); clf; plot(t,theta);
You should see "high" behaviour, so execute
bshigh
1.9749999999999999e+00
1.9937499999999999e+00
1.9843750000000000e+00
Continue the search on your own until you have determined the
critical value, \(\omega_0^\star\)
to at least five decimal places. You can use the function bsfraction
to tell you when you've done so; it outputs
\[ \frac{\omega_{\rm high}- \omega_{\rm low}}{\omega_{\rm current}} \]
so that when it returns a value of about \(10^{-6}\), you'll have determined \(\omega_0^\star\) to about 5 digits. If you type it now you should see
bsfraction
ans =
0.009448818897638
which is about \(10^{-2}\), so you have about 3 or 4 digits to go.
Note that 5 or 6 digits is the minimum precision that you should work out, and, time permitting, you should see how finely you can tune \(\omega_0\) while still being able to discern a difference in the "high" and "low" plots.
If you continue the search to much higher precision, you may eventually
see a type of behaviour which looks like a mixture of "high" or "low".
At that point you are probably exhausting the ability of the simulation
as we have set it up (in particular, with the value of level that
we are using) to narrow in on \(\omega_0^\star\) any further.
It doesn't make much sense to continue the search past that point.
IMPORTANT
Use copy-and-paste to copy the displayed midpoint value and
the assign in to omega0 for the next iteration. It's much
faster and less error prone than manual typing.
IMPORTANT: IF YOU MAKE A MISTAKE DURING THE SEARCH
Ask for help if you don't understand what you are to do, or if you encounter difficulties with the search process.
NOTE
You need to be reasonably attentive and careful when performing
the search. If you accidentally type the wrong command then
use bsundo to back up the search as described in its help
output:
>> help bsundo
bsundo Restores previous search bracket [PHYS210-bs]
Restores the previous bracket, so if you entered bslow when you
should have used bshigh, type
>> bsundo
>> bshigh
and if you entered bshigh when you should have used bslow, type
>> bsundo
>> bslo
Can be applied repeatedly to return to any bracket generated during
the search.
Input arguments
None
Output arguments
None, but prints complete history of search using bsdisp
and then the midpoint value.
Here's an example showing how bsundo works:
>> bsnew(0, 1)
0.0000000000000000e+00
1.0000000000000000e+00
5.0000000000000000e-01
>> bslo
5.0000000000000000e-01
1.0000000000000000e+00
7.5000000000000000e-01
>> bslo
7.5000000000000000e-01
1.0000000000000000e+00
8.7500000000000000e-01
>> bshi
7.5000000000000000e-01
8.7500000000000000e-01
8.1250000000000000e-01
>> bslo
8.1250000000000000e-01
8.7500000000000000e-01
8.4375000000000000e-01
>> bsundo
::::::::::::::
Low values
::::::::::::::
0.0000000000000000e+00
5.0000000000000000e-01
7.5000000000000000e-01
::::::::::::::
High values
::::::::::::::
1.0000000000000000e+00
8.7500000000000000e-01
8.1250000000000000e-01
>> bsundo
::::::::::::::
Low values
::::::::::::::
0.0000000000000000e+00
5.0000000000000000e-01
7.5000000000000000e-01
::::::::::::::
High values
::::::::::::::
1.0000000000000000e+00
8.7500000000000000e-01
>> bsundo
::::::::::::::
Low values
::::::::::::::
0.0000000000000000e+00
5.0000000000000000e-01
::::::::::::::
High values
::::::::::::::
1.0000000000000000e+00
7.5000000000000000e-01
>> bsundo
::::::::::::::
Low values
::::::::::::::
0.0000000000000000e+00
::::::::::::::
High values
::::::::::::::
1.0000000000000000e+00
5.0000000000000000e-01
>> bsundo
Nothing to undo.
Sometimes you may not notice that you have made a mistake until you've performed a few more iterations of the search. This case is often signalled by many consecutive "high" behaviours, or many consecutive "low"s (there's essentially a 50-50 chance of either at any iteration so, for example, if you see 8 "low" values in a row, the chance of that happening legitimately is less than half a percent).
In this case you can use bsundo repeatedly. Each time you
execute it, it will display the low and high values
of the restored bracket followed by the midpoint.
Run
the simulation with each of the displayed low and high values. If you
find that they are correct, i.e. the low value really does
result in "low" behaviour, and the high value yields "high" behaviour,
then continue with the search. Otherwise execute bsundo again
and repeat the process until you have a good bracket.
In the worst case, or if this is all too confusing, you can always start over using
>> bsnew(1.9, 2.05);
which isn't a disaster! It only takes about 50 bisections to narrow the interval to machine precision, a bit over 20 to get it to the minimum precision I've suggested.
Record your estimated value of \(\omega_0^\star\) in
/phys210/$LOGNAME/matlab/README.pendulum
Also, provide answers to the following in that file:
What do you notice about the period of oscillation (for "low" values) as \(\omega_0 \to \omega_0^\star\)
Can you explain what is happening physically?
Make a phase space plot for the "low" value that is closest to your computed critical value using, for example,
plot(theta,omega,'ro');
so that the plot uses small circles rather than a line; recall
that t parameterizes the plot, i.e each circle represents an
instant in time. Save a JPEG version of this plot as
pendulum-phase.jpg in your Matlab directory.
Answer the following in the README.pendulum file.
levelI suggest that you do this only if you have ample free time, which I suspect is unlikely.
Once you've finished with the basic exercise, it would be a much better idea for you to continue with Homework 3!
Repeat the above search procedure for levels 8 and 9 and record
your estimates for
\(\omega_0^\star\)
as a function of level
in your README.pendulum file.
Keep in mind that the specifying level is merely a convenient way
of specifying the "grid spacing", \(\Delta t\),
i.e. the fundamental parameter that governs the accuracy of the FDA.
Note that for each level you will need to reinitialize the search using
bsnew(1.9, 2.05);
Can you make a conjecture about the continuum value of
\(\omega_0^\star\),
i.e. the value we would find for level \(\to\infty\), or
equivalently \(\Delta t \to 0\), and if so, can you prove
that this is what one should expect? Answer in README.pendulum.
Most of you were able to figure out what happens to the pendulum as the critical point is approached. By giving it just the right initial "kick" (value of \(\omega_0\)), we can make it swing around to \(\theta=\pi\) radians (or \(180^\circ\)) where it "balances" in an inverted position relative to the equilibrium position where \(\theta=0\). As we discussed, this is in fact another equilibrium position but of course it is unstable. If the pendulum is perturbed even the slightest amount when it is in that position it will start to swing back towards \(\theta=0\). Nonetheless, it is quite remarkable we can investigate such an exquisitely sensitive phenomenon with such a simple simulation. You will have noticed that as you tuned to \(\omega_0 = \omega_0^\star\) the plots of \(\theta(t)\) developed regions where the curve was flatter and flatter for longer and longer periods of time. This the pendulum "lingering" very near the unstable equilibrium for longer and longer. If we were able to tune infinitely closely to \(\omega_0^\star\), which we can't do because the arithmetic we can do in any computer language is finite precision, we could, in principle make the pendulum balance "upside down" for an infinite amount of time.
This "lingering" is also responsible for the appearance of the phase space plot that some of you were able to make. Because the plot is parametric in time, and because there is a constant amount of time between successive symbols, we see the markers "pile up" at \((\theta=\pi,\,\omega=0)\) and \((\theta=-\pi,\,\omega=0)\).
The optional take-home exercise from the last lab asked you
to investigate the dependence of the critical value of
the initial angular speed of the pendulum, \(\omega_0^\star\),
as a function of the discrete time spacing, \(\Delta t\)
(also known as the time step or
grid spacing or mesh spacing), or equivalently, the level parameter.
If you didn't get a chance to do that, you'll get one in a few minutes, with the scripts that I mentioned above (again please don't run them until I tell you though).
Recall that you did the search in the last lab with level = 7, which
means that there were 65 time steps in the simulation. I asked
you to repeat the search for level = 8 and level = 9 and I
have done that, as well as level = 10, with the following
results:
level omega0* d(omega0*) 4^(level - 7) d(omega0*)
7 1.98348777361164 0.0124290 0.0123
8 1.99591675541708 0.0030651 0.0123
9 1.99898186533854 7.63766e-04 0.0122
10 1.99974563183654
where d(omega0*) = omega0*(level+1) - omega0*(level)
There are two things that are apparent from this table
\(\omega_0^\star\) is converging to some value, and, at the same rate, namely \(O(\Delta t^2)\), as the fundamental dependent variable in the simulation, \(\theta(t^n)\).
Especially given how sensitive the simulation is to changes in the initial conditions as \(\omega_0 \to \omega_0^\star\) this is truly remarkable.
It seems natural to conjecture that the continuum ("exact") value of \(\omega_0^\star\) is exactly 2.
Indeed, from our lecture discussion of the various energy quantities associated with the pendulum we have
\[ T = \frac{\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_{\rm initial} = T_{\rm initial} + V_{\rm initial}
= \frac{(\omega_0^\star)^2}{2} + 0
= \frac{(\omega_0^\star)^2}{2}
\]
\[
E_{\rm final} = T_{\rm final} + V_{\rm final}
= 0 + (1 - \cos(\pi)) = 2
\]
Now, since \(E_{\rm initial} = E_{\rm final}\) we have
\[ \frac{(\omega_0^\star)^2}{2} = 2 \] or \[ \omega_0^\star = 2 \]
This is a nice illustration of something that sometimes happens in computational science (although perhaps not often enough): we tackle a problem which, for some reason or another, we aren't able to get much insight into through "traditional" pen and paper techniques. So we turn to simulation and from the results of the direct numerical calculations we find something "simple" that begs understanding in a more fundamental way. In this case, the appearance of the number 2 in the results is a clear signal that there is probably an "easier" way of getting that answer than through the detailed computations that we made.
In fact, in many areas, including condensed matter physics and astrophysics and applied mathematics, "hard-core theorists" and "numericists" often team-up symbiotically (nonlinearly!). Typically, the numericist will provide results that often are analogous to experimental findings, i.e. something that begs to explained by the theorist. The theorist, gaining some insight from the computational work, can then frequently make predictions of other things that might be found, thus suggesting fruitful avenues for the numericist to explore.
pendulum_searchpendulum_auto_searchpendulum_asxpendulum to compute energies: pendulumeI've coded a modified version of pendulum called pendulume
that computes and returns the energy quantities associated with
the pendulum's dynamics.
This new function has three additional output arguments relative
to pendulum, all of which are row vectors of length nt (the
same length as omega and theta):
T -> Kinetic energy as a function of discrete time
V -> Potential energy as a function of discrete time
E -> Total energy as a function of discrete time
As usual, we can view the definition of the function in the command
window using the type command
>> type pendulume
function [t, theta, omega, T, V, E] = pendulume(tmax, level, theta0, omega0, ...
tracefreq)
% pendulume Solves the nonlinear pendulum equation using an O(deltat^2) FDA [phys210-pendulum]
%
% Variant of pendulum which computes energy quantities
%
% Input arguments
%
% tmax: (real scalar) Final solution time.
% level: (integer scalar) Discretization level.
% theta0: (real scalar) Initial angular displacement of pendulum.
% omega0: (real scalar) Initial angular velocity of pendulum.
% tracefreq: (optional integer scalar) Frequency of tracing output,
% 0 disables tracing.
%
% Output arguments
%
% t: (real vector) Vector of length nt = 2^level + 1 containing
% discrete times (time mesh).
% theta: (real vector) Vector of length nt containing computed
% angular displacement at discrete times t(n).
% omega: (real vector) Vector of length nt containing computed
% angular velocity at discrete times t(n).
% 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).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tracing control: if 5th arg is supplied base tracing on that input,
% otherwise use local defaults.
if nargin > 4
if tracefreq == 0
trace = 0;
else
trace = 1;
end
else
trace = 1;
tracefreq = 100;
end
if trace
fprintf('In pendulume: Argument dump follows\n');
tmax, level, theta0, omega0
end
% Define number of time steps and create t, theta and omega arrays of
% appropriate size for efficiency.
nt = 2^level + 1;
t = linspace(0.0, tmax, nt);
theta = zeros(1, nt);
omega = zeros(1, nt);
% Determine discrete time step from t array.
deltat = t(2) - t(1);
% Initialize first two values of the pendulum's angular displacement.
theta(1) = theta0;
theta(2) = theta0 + deltat * omega0 - 0.5 * deltat^2 * theta0;
if trace
fprintf('deltat=%g theta(1)=%g theta(2)=%g\n',...
deltat, theta(1), theta(2));
end
% Initialize first value of the angular velocity.
omega(1) = omega0;
% Evolve the oscillator to the final time using the discrete equations
% of motion. Also compute an estimate of the angular velocity at
% each time step.
for n = 2 : nt - 1
% This generates tracing output every 'tracefreq' steps.
if rem(n, tracefreq) == 0
fprintf('pendulume: Step %g of %g\n', n, nt);
end
theta(n+1) = 2 * theta(n) - theta(n-1) - deltat^2 * sin(theta(n));
omega(n) = (theta(n+1) - theta(n-1)) / (2 * deltat);
end
% Use linear extrapolation to determine the value of omega at the
% final time step.
omega(nt) = 2 * omega(nt-1) - omega(nt-2);
% Compute vectors of kinetic, potential and total energy.
% Note how straightforward the modifications of pendulum are due
% to octave's ability to perform whole-array operations, but also
% observe that we need to use '.*' in the computation of T to ensure
% that the multiplication is done element-by-element.
T = 0.5 * omega .* omega;
V = 1 - cos(theta);
E = T + V;
end
pendulume: tpetpe is a driver script for pendulume that does the following
Performs a basic test of pendulume and plots all three energy
quantities.
Performs a convergence test using level = 6, 7, 8, 9 that plots
deviations in the total energy relative to the initial time,
\(t=0\).
Performs another convergence test using the same levels but where the plotted deviations have been scaled so that near-coincidence of the graphs signals the expected \(O(\Delta t^2)\) convergence to conservation.
Let's give it a whirl ...
>> tpe
The script demonstrates that the simulations (calculations) do
exhibit energy conservation: not exactly, but in precisely the same sense
that the numerical values
\(\theta^n \equiv\theta(t^n)\) that pendulum and
pendulume compute satisfy the equations of motion for the
pendulum: i.e. the errors go to 0 like \(\Delta t^2\).
There's an important point to make here.
As physicists, we have the fundamental importance of conservation of energy (and other conserved quantities like momentum and angular momentum) drilled into us from the moment we start studying physics seriously.
So is natural for a physicist doing a simulation like this to focus on the issue of energy conservation and perhaps worry if he/she sees that, as is the case here, it's not exactly conserved.
However, the important point is that the energies that we are computing are derived values; we calculate them "after the fact" once we have determined the quantity that is truly of fundamental importance: namely the dynamical variable \(\theta^n\equiv\theta(t^n)\). Once we have evidence that \(\theta^n\) is converging at the expected rate, then if we didn't observe convergence in energy conservation as well, it would tell us not that energy really wasn't being conserved, but that we had computed one or more of the energies incorrectly (e.g. we had used an incorrect formula to compute one or more of \(T\), \(V\) and \(E\). In fact, some of you may encounter this phenomenon when you work the exercise coming up.
If the solution itself is converging then there is simply no way that the the total energy cannot converge to conservation.
In other words, it's fine to monitor conserved quantities, but it's much more important to study the convergence behaviour of the dynamical variables themselves.
lpendulum to compute energies: lpendulume.mCopy ~phys210/matlab/lpendulum.m to your Matlab
directory /phys210/$LOGNAME/matlab and rename the copy
lpendulume.m.
Modify lpendulume.m so that it defines a function, lpendulume, with
a header
function [t theta omega T V E] = lpendulume(tmax, level, theta0, omega0)
that does precisely what lpendulum.m does, but that also includes a
calculation of the vectors
T, V and E containing the kinetic,
potential and total energies, respectively, for the simulation. Each
of the vectors should be the same length as t, theta and omega,
i.e. they should each have nt elements, where nt is defined in
terms of level in the usual way (you won't have to change the
definition).
In short, lpendulume.m is to lpendulum.m as pendulume.m is to pendulum.m
IMPORTANT
Be sure to use the correct expressions (formulae) for the kinetic and especially the potential energy. The latter is not the same as it is for the nonlinear pendulum.
If needed, refer to page 32 of the pendulum notes via the syllabus section of the course page) for the relevant formulae for \(T\), \(V\) and \(E\).
tpe to test lpendulume: tlpe.mNow, copy ~phys210/matlab/tpe.m to your Matlab
directory and rename it tlpe.m.
Note that tpe.m is a script, i.e. it has no function header, and
tlpe.m will be a script as well.
Modify tlpe.m so that it invokes lpendulume rather than pendulume,
and then execute the script to demonstrate convergence of the
energy quantities for the linear simulation
Do you notice anything special/peculiar about the behaviour of the convergence plots (and, in fact, all of the plots) at late times (i.e. at the very end of the simulations) and, if so, can you hypothesize what might be causing the behaviour?
Add your observations/answers to the README.pendulum file you
created in the last lab.