------------------------------------------------------------------------- PHYS210 NONLINEAR PENDULUM LAB 2 SUMMARY 1) Discussion of the critical value, omega0*, for the nonlinear pendulum 2) Modification of pendulum.m (-> pendulume.m) to compute energy quantities 3) Convergence test of total energy conservation of nonlinear pendulum 4) Exercise: modification of lpendulum.m (-> lpendulume.m) to compute energy quantities, modification of tpe.m -> tlpe.m to convergence test energy conservation 5) 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 - 7) d(omega0*) 7 1.983487769961357 0.0124290 0.0124 8 1.995916742086411 0.0030651 0.0123 9 1.998981863260269 7.63768e-04 0.0122 10 1.999745631217956 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 our 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 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Start Matlab - Look at the definition of pendulume.m >> type pendulume function [t theta omega T V E] = pendulume(tmax, level, theta0, omega0) % pendulume Solves nonlinear pendulum equation using second order FDA % % 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. % % 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). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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 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 %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); % 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 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3) Use driver script 'tpe.m' to illustrate computation of energy quantities and their convergence behaviour +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ >> tpe +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 4) EXERCISE: +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ------ Part 1 ------ - Copy ~phys210/matlab/lpendulum.m to ~/matlab, renaming the file as lpendulume.m in the process, i.e., in a terminal window execute % cd ~/matlab % cp ~phys210/matlab/lpendulum.m 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 - I.e. lpendulume.m is to lpendulum.m as pendulume.m is to pendulum.m - If needed, refer to page 32 of the pendulum class notes (available via the syllabus section of the course page) for the relevant formulae for T, V and E ------ Part 2 ------ - Copy ~phys210/matlab/tpe.m to ~/matlab, renaming the file as tlpe.m in the process, i.e., in a terminal window execute % cd ~/matlab % cp ~phys210/matlab/tpe.m tlpe.m - Note that tlpe.m is a SCRIPT file, i.e. it has no function header - 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? +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 5) Free time to work on Homework 3 and projects! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++