------------------------------------------------------------------------- PHYS210 NONLINEAR PENDULUM LAB 1 SUMMARY NOTE: Item 1) below can be worked "automatically" via the instructor supplied script, 'tp'. I.e. type 'tp' at the Matlab prompt. 1) Studying / using instructor-supplied Matlab 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 7) 4) EXERCISE: Investigating dependence of omega0* on level ------------------------------------------------------------------------- ------------------------------------------------------------------------- - Solution to last lab's exercise available in ~phys210/matlab/demoreturn.m - Start Matlab: For the lab exercises, it should now not matter from where you start Matlab, including the launcher. However, any .m files that you create during the labs should be located in ~/matlab. This will ensure that Matlab can always find them when you type the corresponding function/script name - If execution of any of the functions or scripts below does not work, it probably means that your Matlab path, that we discussed and modified in the last lab session, is not set properly. In any case, if you encounter any problems, ask for help immediately! - ADDITIONAL HINT for HW3, PROBLEM 2: demonstration function, 'mytranspose', that shows use of nested for loops for defining elements of a 2D array (matrix), as well as the mechanism for determining the dimensions of an array - Display the definition of 'mytranspose' >> type mytranspose function mattrans = mytranspose(mat) % Transposes a matrix as a demonstration of the use of nested for loops % for array operations. % % Input argument: mat --- m x n matrix (2-D array) % % Output argument: mattrans ---- n x m matrix, transpose of mat % Determine dimensions of array: recall that size(a) returns % a row vector of the dimensions. sizemat = size(mat); % First element is the number of rows. m = sizemat(1); % Second element is the number of columns. n = sizemat(2); % Create the output array, initialized to zero. % Preallocating an array in this way is more efficient % than "growing" the array element by element. mattrans = zeros(n, m); % Loop over all elements of mat using nested for's for ii = 1 : m for jj = 1 : n % Transpose individual elements ... mattrans(jj, ii) = mat(ii, jj); end end end - Test the function >> m = [11:13; 21:23; 31:33; 41:43] m = 11 12 13 21 22 23 31 32 33 41 42 43 >> mytranspose(m) ans = 11 21 31 41 12 22 32 42 13 23 33 43 - You should also note the additional point concerning the use of the ';' as a statement terminator to suppress output - Specifically, observe that IT IS NOT NECESSARY TO use semi-colons at the end of lines that constitute a PORTION of a control structure. - For example, we write, for example ... if a > b val = a + b; else val = b - a; end ... NOT ... if a > b; val = a + b; else; val = b - a; end; - Similarly, use, e.g. ... for i = 1 : 10 vec(i) = i^2; end ... NOT ... for i = 1 : 10; vec(i) = i^2; end; - Although the additional semi-colons in the above code snippets ARE valid syntactically, they are superfluous ------------------------------------------------------------------------- +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1) Studying / using Matlab functions to solve nonlinear & linear pendulum equations +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - From within the Matlab command window, examine the definition of 'pendulum' which solves the nonlinear pendulum equation using O(h^2) finite difference techniques as discussed in class. >> 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 - There is a precisely analogous function, lpendulum, which solves the *linear* pendulum equation (i.e. the one that you study in elementary physics courses). >> 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 - Now call pendulum, being careful to specify a VECTOR of three names on the right 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 - Call -> lpendulum <- with the same arguments, but being sure to use DIFFERENT NAMES ON THE LEFT HAND SIDE FOR THE ANGULAR DISPLACEMENT AND VELOCITY (ltheta and lomega) >> [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' --- 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---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 >> 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 -------------------- 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 tmax, theta0 and omega0 are held fixed, while the grid spacing is varied - For this test we will use levels 6, 7 and 8 - FROM THIS POINT ONWARD, I RECOMMEND THAT YOU USE "copy and paste" TO ENTER THE COMMANDS, AND TO EXPEDITE THIS I WILL NOT PUT THE PROMPT IN FRONT OF EACH MATLAB COMMAND LINE, BUT RATHER ON LINES BEFORE AND AFTER EACH COMMAND OR BLOCK OF COMMANDS - 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 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 >> theta7 = theta7(1:2:length(theta7)); theta8 = theta8(1:4:length(theta8)); >> - 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) a direct measure of the error in the numerical solution - Now scale (multiply) dtheta78 by 4 and replot 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, deltat---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 --- 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 delt -> 0 (i.e. as level -> infinity), the leading order---O(delt^2)---error term in the finite difference solution will become increasingly dominant over the higher order terms ( O(delt^4), ...) - We could continue the convergence test to even higher levels, level = 10, 11, ... , but we should note that, for sufficiently high level (sufficiently small delt), 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 is plotted against theta (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 - 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 >> 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('"lo(w)" and "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 a.k.a. BISECTION. - ASIDE: As remarked 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 look at the contents of 'tp.m', which is located at ~phys210/matlab/tp.m as an example of a basic script (use 'type tp' in Matlab, or 'more ~phys210/matlab/tp.m' in a terminal window). ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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 you are still confused about this topic after our discussion of it They are 1) bsnew 2) bscurr 3) bslo 4) bshi 5) bslohi 6) bsfrac With usages and functions as follows 1) bsnew usage: bsnew FUNCTION: Initializes a new binary search in the current directory with initial interval [, ] 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 a terminal window open, and ensure that the working directory is ~/matlab % cd ~/matlab - As you do this numerical experiment---a search for the critical value omega0*---you will be executing commands in BOTH the terminal window (i.e. at the shell prompt, denoted '%' as usual) and the Matlab command window (i.e. at the Matlab prompt, denoted '>>' as usual). Be careful that you type each command in the proper window - Initialize the search: Note that the script returns the values defining the initial interval - In the terminal window ... % bsnew 1.9 2.05 1.9 2.05 - NOTE: If you see the error message -bash: bsnew: command not found your bash path is not properly set---ask for help! - Get the current midpoint of the bracket (i.e. the average of the "lo" and "hi" values) % bscurr 1.975000000000000e+00 - In the Matlab window ... - Change the format to long so that we can see values to full precision >> format long >> - Set the fixed problem parameters - NOTE: WE WILL PERFORM THE SEARCH USING level 7 >> tmax = 40; level = 7; theta0 = 0; >> - Set omega0 to the midpoint of the interval, i.e. to the value returned by the script 'bscurr' >> 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 oscillating (i.e. "lo" behaviour), so in THE TERMINAL WINDOW EXECUTE ... % bslo 1.975000000000000e+00 2.05 - Note how the updated interval values are echoed. - Now, generate the new midpoint value ... % bscurr 2.012500000000000e+00 - Return to the Matlab window, and repeat the simulation with the new value of omega0 >> omega0 = 2.012500000000000e+00 [t theta omega] = pendulum(tmax, level, theta0, omega0); clf; plot(t,theta); >> - This time you should observe "hi" behaviour---i.e. theta growing with time, so return to the shell window and execute ... % bshi 1.975000000000000e+00 2.012500000000000e+00 - And then generate the new midpoint value .. % bscurr 1.993750000000000e+00 - Go back to Matlab, set omega0 to this value, and run the simulation >> omega0 = 1.993750000000000e+00 [t theta omega] = pendulum(tmax, level, theta0, omega0); clf; plot(t,theta); >> - You should see hi behaviour - In the shell window, execute ... % bshi 1.975000000000000e+00 1.993750000000000e+00 % bscurr 1.984375000000000e+00 ------------------------------------------------------------------------ EXERCISE ------------------------------------------------------------------------ - CONTINUE THE SEARCH ON YOUR OWN UNTIL YOU HAVE DETERMINED THE CRITICAL VALUE, omega0*, TO AT LEAST FIVE DECIMAL PLACES, I.E. SO THAT WHEN YOU EXECUTE ... % bsfrac ... at the shell prompt, you get a number which is APPROXIMATELY 1e-6 - ASK FOR HELP IF YOU DON'T UNDERSTAND WHAT YOU ARE TO DO, OR IF YOU ENCOUNTER DIFFICULTIES WITH THE SEARCH PROCESS - Record your computed value of omega0 in the file ~/matlab/README.pendulum - NOTE: BE VERY CAREFUL DOING THE SEARCH. IT IS EASY TO "LOSE THE BRACKET" SO THAT EITHER OF THE CURRENT "HI" OR "LO" VALUES ARE NOT REALLY "HI"/"LO" ... IF THIS APPEARS TO HAVE HAPPENED TO YOU ASK FOR HELP --- ALTERNATIVELY, USE THE ARROW KEYS TO RECALL THE VALUES OF omega0 YOU HAVE USED, AND ATTEMPT TO FIND THE MOST RECENTLY USED "HI" AND "LO" VALUES (I.E. VALUES THAT *DO* SHOW THE "HI" AND "LO" BEHAVIOUR), THEN REINITIALIZE THE SEARCH BY EXECUTING % bsnew IN THE TERMINAL WINDOW ------------------------------------------------------------------------ QUESTIONS ------------------------------------------------------------------------ - WHAT DO YOU NOTICE ABOUT THE PERIOD OF OSCILLATION (for "lo" values) AS omega0 -> omega0*? - CAN YOU EXPLAIN WHAT IS HAPPENING PHYSICALLY? - MAKE A PHASE SPACE PLOT FOR THE "lo" VALUE THAT IS CLOSEST TO YOUR COMPUTED CRITICAL VALUE. USE >> plot(theta,omega,'ro'); >> so that the plot uses small circles rather than a line; recall that t parametrizes the plot, i.e each circle represents an instant in time - WHAT DO YOU OBSERVE, AND CAN YOU EXPLAIN IT PHYSICALLY? +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 4) Investigating dependence of omega0* on level ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ------------------------------------------------------------------------ EXERCISE ------------------------------------------------------------------------ - Repeat the above search procedure for levels 8 and 9 and record your estimates for omega0* as a function of level in the file README.pendulum in ~/matlab. - 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 omega0* (i.e. THE VALUE FOR level -> infinity, deltat -> 0) AND, IF SO, CAN YOU PROVE THAT THIS IS WHAT ONE SHOULD EXPECT??