Physics 6730 Assignment 3

For each of the following exercises create the specified files with your answer(s). Submit your homework using the course submit utility.

Exercise 1. Anharmonic oscillator: Simple 4th order Runge-Kutta

This is a programming exercise to solve the differential equation for a damped, forced true pendulum (anharmonic oscillator). This system is particularly interesting, since it exhibits both periodic and chaotic behavior, depending on the choice of parameters. See the handout Chaos for a brief discussion of the problem. You are to find the solution to the pair of equations:

     dy1/dt = y2
     dy2/dt = -sin(y1) - y2/q + g cos(wd t)

The solution to the above equations is not known exactly, so it must be found numerically. Part of the challenge is to test for convergence. We will be applying three numerical techniques to solve the problem, so we can compare their accuracy against the known exact solution. The methods we will be comparing are these:

1. Runge Kutta 4th order
2. Runge Kutta 4th order with adaptive step-size control
3. Modified midpoint with rational function extrapolation.

In this exercise you will write ONLY the routine that evaluates the rhs of the differential equation AND the Makefile. The programs you will bundle to hand in are these:

File	Purpose of procedure
anhrk4.cc	Driver (main) program. Calls NR::rk4, anhderivs.
rk4.cpp	Does one Runge-Kutta 4th order step. Note that all Numerical Recipes routines are placed in the NR namespace, so the full entry point name of the subroutine is NR::rk4.
anhderivs.cc	Evaluates the rhs of the ODEs.
Makefile	Rule for compiling the above.

The file anhrk4.cc is provided in ~p6730/examples/a03. You should copy it. The file rk4.cpp is in ~p6730/examples/nmr. You don't need to copy it. The object file rk4.o is in ~p6730/lib/libnmr.a and is linked when you use the linker flags -L${INSTRUCTOR}/lib -lnmr. The environment variable ${INSTRUCTOR} is set to the instructor's home directory.

In summary, you are to write only one of the above procedures and the make file. Most of the work in this exercise is in providing the proper interface with the stepping routines, and in understanding what the routines are doing.

All computations are to be done in DOUBLE precision.

1. The subprogram file rk4.cpp.

   Please see Numerical Recipes for details. The source code is in ~p6730/examples/nmr. Note that rk4 is a simple stepper. It expects that y, dydt, and t are set to the correct values at the beginning of the interval. It returns the new values in yout, computed using the 4th order Runge Kutta method. It DOES NOT step the variable t, DOES NOT call anhderivs to get dydt at the new values, and DOES NOT actually replace the old values y with the new values. That must be done in the main program.
   
   

2. The main program file anhrk4.cc

   Note that all interactive output is to standard error. Standard output has a list of values

        t   theta(t)   omega(t)
   

   
   
   
3. Your subprogram file anhderivs.cc

   This subprogram is called by rk4 and by anhrk4. The purpose of this routine is to take as input the values y(1) and y(2) and t and evaluate the rhs of the ODEs, returning the result in dydt(1) and dydt(2). It may be coded in just a few lines.

   The values of the parameters and the number of calls are kept in the params class object p. See the main program anhrk4.cc for the definition of the params class. This variable is defined as a global. Storage is allocated in the main program, so you should declare it as extern params p in your subprobram.
   
   

4. Your Makefile as usual.

Exercise 2. Running the code

Note: for small dt, the volume of output grows considerably. You can select the very last values by piping stdout into tail.

For amusement, you might try making a phase space plot of the output by selecting the pair (theta, omega):

   anhrk4 < data > plotfile

Please put your analysis of the global error in a file rk4.txt

Exercise 3. Quality control

Here are the changes:

File	Purpose of procedure
anhrkqc.cc	Driver (main) program. Calls NR::rkqs, anhderivs.
rkqs.cpp	Does one RKQC step
rk4.cpp	Does one Runge-Kutta 4th order step
anhderivs.cc	Evaluates the rhs of the ODEs.
Makefile	Rule for compiling the above.

In summary, hand in anhrkqc.cc and Makefile. Your main program anhrkqc.cc and Makefile are just modifed versions of the ones from Exercise 1. The Makefile you submit should compile the targets anhrk4 for Exercise 1 and anhrkqc for this exercise. The file anhderivs.cc should work for both of them.

As before, all computations are to be done in DOUBLE precision.

1. The supplied program NR::rkqs

   See the comments in the code and Numerical Recipes for a definition of the variables and a detailed discussion of rkqs.cpp. To summarize, the way it works is that you tell it you want to step a distance dt. It tries to do that and it estimates the presulting error in y(1) and y(2). It compares the error with the maximum allowed errors that you specify for each, namely, yscal(1)*eps and yscal(2)*eps. If either estimated error is larger than allowed, it tries again with a shorter step. In that case the output value dtdid will be less than dt. If the estimated errors are both smaller than allowed, it figures that you could get by with a larger step. In that case it takes the step dt that you asked for, but suggests a longer step dtnext for the next time. You are free to accept its advice and set dt equal to dtnext, or you may give it a different value dt.

   The program expects that y, dydt, and t are set to the correct values at the beginning of the interval. It returns the new values in y, computed using the adaptive 4th order Runge Kutta method, and steps the variable t, but DOES NOT call anhderivs for dydt at the new values. Thus your main program must call anhderivs each time it calls rkqs. Unlike rk4.cpp, however, it DOES update the independent variable t and it DOES update y.
   
   

2. Modifications to the main loop in anhrk4 to create anhrkqc
   • You will have to read one more parameter, eps.
   • yscal(i)*eps determines the maximum allowed error in y(i) The following choice seems to work well for this problem:

              yscal(1) = fabs(y(1)) + fabs(dydt(1)*dt) + 1.e-8
              yscal(2) = fabs(y(2)) + fabs(dydt(2)*dt) + 1.e-8
     

     The array yscal should be recomputed BEFORE each call to NR::rkqs.
   • After each call to NR::rkqs, replace dt with dtnext. That is, accept the suggestion.
   • You should design the loop over calls to NR::rkqs so that you quit the loop if the next step would give t > tmax. However, since we want to integrate up to t = tmax exactly, you will then have to make one final call to NR::rkqs with dt = tmax - t to complete the solution.

Exercise 4. Pick the winner!

Notice that with anhrk4 you have only the step size dt to play with. With the other methods, you adjust eps to control the error. As you decrease these tolerances, the method has to call anhderivs more times and the answer becomes more accurate until roundoff errors take over.

Run all three programs anhrk4, anhrkqc and anhbust with the same input parameters as in Exercise 2, and adjust dt and eps (as appropriate) so that the error in theta at t = 200 is less than .001.

For each of the three methods vary the tolerance until you find the largest value of either eps or dt that gets the correct answer to four digits. (Rounding is permissible.) To shorten your work, vary eps only by factors of 10 and vary dt in the sequence .03, .01, .003, etc. In your answer file make a table giving the score (number of calls to anhderivs) for each method. Also give the value of theta and omega at t = 200.

Please put your results and answers a file rkqc.txt.