-------------------------------------------------------------------------
PHYS210 2009-11-03 LAB SUMMARY

NOTE: Item 1) below can be worked "automatically" via the instructor
supplied script, 'tp'.  I.e. type 'tp' at the octave prompt.

1) Studying / using instructor-supplied octave functions to solve
   nonlinear & linear pendulum equations

2) Use of instructor-supplied bash scripts to expedite binary search

3) Finding critical omega0 value, omega0*, for nonlinear pendulum using 
   binary search (level 6)

4) Investigating dependence of omega0* on level

NOTE: In the following, '>>' denotes the octave prompt

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1) Studying / using instructor-supplied octave functions to solve
   nonlinear & linear pendulum equations
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Change to your ~/octave directory and start octave

% cd ~/octave 
% octave

>> % Examine definition of 'pendulum' which solves the nonlinear 
>> % pendulum equation using O(h^2) finite difference techniques 
>> % as discussed in class.  Recall that 'type <function name>'
>> % will display the definition of a function, as well as the 
>> % location of the .m file

>> % Disable pager
>> more off

>> type pendulum

pendulum is the user-defined function defined from: /home/phys210/octave/pendulum.m

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.
%
%  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).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % 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
      pause(1)
   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);

   % Intialize 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

>> % There is a similar function, lpendulum, which solves the 
>> % *linear* pendulum equation (i.e. the one that you study in 
>> % elementary physics courses.

>> type lpendulum

lpendulum is the user-defined function defined from: /home/phys210/octave/lpendulum.m

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.
%
%  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).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % 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
      pause(1)
   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);

   % Intialize 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

      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

>> % Try using the functions:  in all cases below, we will start with
>> % theta0 = 0 (i.e. so that the pendulum is hanging vertically), and
>> % will vary the initial angular velocity omega0

>> % Start with a small value of omega0, and compute with input arguments
>> % as follows

>> tmax = 40.0; theta0 = 0; level = 8; omega0 = 0.01;

>> [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

>> [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
>> clf; hold on; plot(t, theta, 'r'); plot(t, ltheta, '.og');

>> % Observe how there is very little difference between the results, 
>> % as one would expect.

>> % Perform a basic convergence test for the nonlinear pendulum using
>> % levels 6, 7 and 8

>> [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 from the 3 calculations on the same plot
>> clf; hold on; plot(t6, theta6, 'r-.o'); 
>> plot(t7, theta7, 'g-.+'); plot(t8, theta8, 'b-.*');

>> % "Filter" 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)

>> theta7 = theta7(1:2:length(theta7));
>> theta8 = theta8(1:4:length(theta8));

>> % Compute the differences in the grid functions from level to level ...

>> dtheta67 = theta6 - theta7;
>> dtheta78 = theta7 - theta8;

>> % ... and make a plot of those functions

>> clf; hold on; plot(t6, dtheta67, 'r-.o'); plot(t6, dtheta78, 'g-.+');

>> % Now scale dtheta78 by 4 and replot with dtheta67.  
>> % If the convergence is second order, the curves should be 
>> % nearly coincident

>> 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

>> % Now use a sequence of ever larger values of omega0 at level 8
>> % Remember, this means that we will be giving the pendula ever 
>> % increasing initial angular velocities, so the resulting amplitude
>> % of the oscillation will also increase.
>> %
>> % A key result that you can see immediately is whereas the period 
>> % (frequency) of the linear pendulum (green) stays constant (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)
>>
>> tmax = 40; level = 8; theta0 = 0;
>> 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);
	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');
>> % Can you explain this result physically?

>> % 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('"lo(w)" \& "hi(gh)" 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, or bisection.

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2) Use of instructor-supplied bash scripts to expedite binary search
   (also known as the "bisection method"
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
There are a set of scripts in ~/phys210/bin (which should therefore 
be in your path) that can be used to expedite a binary search (google
"bisection method" if this concept is unfamiliar to you)

They are 

1) bsnew 
2) bscurr
3) bslo
4) bshi
5) bslohi
6) bsfrac

With usages and functions as follows

1) bsnew 

usage: bsnew <lo> <hi>

FUNCTION: Initializes a new binary search in the current directory 
with initial interval [<lo>, <hi>]

2) bscurr

usage: bscurr

FUNCTION: Returns current midpoint of interval

3) bslo

usage: bslo

FUNCTION: Updates interval so that current midpoint becomes new "lo" value

4) bshi

usage: bshi

FUNCTION: Updates interval so that current midpoint becomes new "hi" value

4) bslohi

usage: bslohi

FUNCTION: Echoes interval end-point values ("lo" and "hi")

5) bsfrac

usage: bsfrac

FUNCTION: Returns normalized width of current interval.  Tells you to 
how many digits you have narrowed the search.  When bsfrac returns 
something of order 10^{-16} you have narrowed the search about as far
as you can

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3) Finding critical omega0 value, omega0*, for nonlinear pendulum using
   binary search (level 6)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Have at least two terminal windows open, and ensure that you are in
your ~/octave directory in both

Use one window to use the above scripts to manage the binary search, 
and the other to run pendulum with the value returned by bscurr for 
omega0 

WINDOW 1

% cd ~/octave
% # Initialize the search: Note that the script returns the 
% # values defining the initial interval
% bsnew 1.9 2.05
1.9
2.05

% # Get the current midpoint
% bscurr
1.975000000000000e+00

WINDOW 2 

% cd ~/octave
% octave

>> format long

>> % Set fixed problem parameters 
>> tmax = 40; level = 6; theta0 = 0;

>> % Set omega0 to midpoint
>> omega0 = 1.975000000000000e+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 growing with time (i.e.
>> % "hi" behaviour, so in WINDOW 1 execute

WINDOW 1

% bshi
1.9
1.975000000000000e+00

% # ... note how the updated interval values are echoed.
% # Generate the new midpoint value ...

% bscurr
1.937500000000000e+00

% # ... and go back to WINDOW 2 (the octave session), repeat the 
% # process with the new value of omega0

WINDOW 2

>> omega0 = 1.937500000000000e+00
>> [t theta omega] = pendulum(tmax, level, theta0, omega0); clf; plot(t,theta);

>> # Observed hi behaviour

WINDOW 1 

% bshi; echo; bscurr
1.9
1.937500000000000e+00

1.918750000000000e+00

WINDOW 2 

>> omega0 = 1.918750000000000e+00;
>> [t theta omega] = pendulum(tmax, level, theta0, omega0); clf; plot(t,theta);

>> # Observed lo behaviour 

WINDOW 1

% bslo; echo; bscurr

1.918750000000000e+00
1.937500000000000e+00

1.928125000000000e+00

AND CONTINUE UNTIL YOU HAVE DETERMINED THE CRITICAL VALUE, omega0*, TO AT LEAST 
SIX DECIMAL PLACES

WHAT DO YOU NOTICE ABOUT THE PERIOD OF OSCILLATION (for "lo" values)
AS omega0 -> omega0*?

CAN YOU EXPLAIN WHAT IS HAPPENING PHYSICALLY?

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4) Investigating dependence of omega0* on level
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Repeat the above search procedure for levels 7 and 8 and record 
your estimates for omega0* as a function of level in README.pendulum
in ~/octave.

CAN YOU MAKE A CONJECTURE ABOUT THE CONTINUUM VALUE OF omega0* (i.e.
THE VALUE FOR level -> infinity, deltat -> 0) AND, IF SO, CAN YOU 
PROVE THAT THIS IS WHAT ONE SHOULD EXPECT??