------------------------------------------------------------------------- PHYS210 NONLINEAR PENDULUM 2 AND N-BODY LAB 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) Demonstration of xfpp3d, utility for visualization of N-body results, and nbodyout function to write data to a file in the format needed by xfpp3d 5) IMPORTANT note concerning use of random number function (rand) for generating initial positions and velocities in N-body codes 6) 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 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 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % cd octave % octave >> type pendulume function [t theta omega T V E] = pendulume(tmax, level, theta0, omega0) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Solves nonlinear 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). % 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); % 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('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) Demonstration of xfpp3d, utility for visualization of N-body results +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - IMPORTANT: xfpp3d must be executed from the shell, NOT in octave - In a terminal window, change to your octave directory % cd % cd octave - Copy the sample data file ~phys210/octave/nbody3.dat % cp ~phys210/octave/nbody3.dat . - Invoke xfpp3d with the '-h' option to get help for the utility (also available via Home Page -> Online Course Resources -> Visualization Utilities -> xfpp3d % xfpp3d -h usage: xfpp3d [-h] [-c] [-s] [-w] Displays animations of three dimensional particle motion. Options: -h Displays this help message. -c Enables user-specified coloring (RGB) of particles. -s Renders particles as solid spheres. -w Renders particles as wireframe spheres. xfpp3d reads ASCII input from standard input in one of the following two formats: (Note that ! denotes the start of a comment.) . . . All values are expected to be real numbers except for which is an integer. - Execute xfpp3d using the sample data file as input, and supplying the option -c so that the particles are coloured IMPORTANT!! At the current time, xfpp3d only takes input from the standard input, so USE INPUT REDIRECTION ('<') % xfpp3d -c < nbody3.dat - I will demonstrate usage of xfpp3d ------------------------------------------------------------------------ - There are functions in ~phys210/octave that can be used to output N-body data to a file in the format that xfpp3d requires nbodyout nbodyout1 - In an octave session execute >> help nbodyout >> help nbodout1 for usage information and >> type nbodyout >> type nbodyout1 to see the function definitions - As with the other instructor-supplied scripts, and assuming that you are working on the lab computers, you should be able to use these directly in your programs +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 5) IMPORTANT note concerning use of random number function (rand) for generating initial positions and velocities in N-body codes +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - First, note that as was discussed in class, and as is described in the N-body notes http://laplace.phas.ubc.ca/210/Doc/fd/nbody.pdf I recommend that you test gravitational N-body codes using a configuration describing two particles in mutual circular orbits about their center of mass. This set-up should allow you to investigate convergence properties of your implementation, including the convergence of conservation of total energy. - If you wish to experiment with random initial positions and velocities, be aware that, under normal circumstances, each time you run your script, a NEW (i.e. distinct) sequence of random numbers will be generated, so that you will NOT be able to perform a convergence test by re-running the calculation with varying values of level (since the initial conditions will NOT be the same from run to run) - If you DO wish to perform convergence tests with random initial values, then you can do so by setting the 'state' of the random number function 'rand' prior to using it to generate the actual random numbers - There is a demonstration script that illustrates this facility ~phys210/octave/trand that you should be able to view using >> type trand and execute via >> trand - Use >> help rand for more information if needed, or ask me for help.