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