PHYCS 6730 Lab Exercise: Ordinary Differential Equations: Global Error

In these exercises we will be comparing the error as a function of step size for three methods: 4th order Runge Kutta (anhrk4), 4th order Runge Kutta with quality control (anhrkqc), and 4th order Runge Kutta with Bulirsch-Stoer extrapolation (anhbust). The methods are described in Numerical Recipes by Press et al.

The executable code you need is in ~p6730/exercises/ode. It is already compiled.

This code solves the differential equation 

     d2y/dt2 = -sin(y) - 1/q dy/dt + g cos(wd t).
The variable y is the angle of the pendulum from vertical (downward) equilibrium. dy/dt is the angular velocity. The coefficients q, g, and wd come respectively from the damping, driving strength (not to be confused with gravitational acceleration), and driving frequency. The initial conditions are specified by the starting angle and angular velocity at time t = 0. We are interested in finding the angle at some later time Tmax. We might also be interested in knowing the angle at all times between t = 0 and t = Tmax.

Exercise 1 Chaotic Solutions

Ordinarily the damped, driven pendulum exhibits a transient oscillation which settles down to steady harmonic motion at the driving frequency. When the pendulum is driven at a suitably large amplitude and with appropriate damping, however, it begins to show erratic behavior. To see such a solution, run anhrk4, redirecting input parameters from the file data_chaos, which you will find in the same directory. The input parameters are, in order, 
  q g wd
  y(0) dy/dt(0)
  Tmax dt
  eps
where eps controls the desired accuracy (used only for anhrkqc and anhbust).

So you can visualize the result, run it like this: 

    anhrk4 < data_chaos > out.rk4
The output file out.rk4 has three values per line. The first is the current time, the second, the current value of y, and the third, the current value of dy/dt. Use gnuplot to plot y vs t.

Describe what you see. How does this behavior differ from the usual transient decay to a steady state behavior for a driven, damped oscillator?

Exercise 2. Effect of quality control

Next, try running anhrkqc and anhbust with the same input data. Give your output files distinctive names: out.rkqc and out.bust.

What differences do you see between the results? What are the final values of y and dy/dt? Report these numbers in your answer file. How many steps were taken? Report these numbers. What about the number of calls to anhderivs, which evaluates the rhs of the ODE? This is a measure of the computational effort. Report this number. Try comparing two plots by using the gnuplot command 

 plot
'out.bust' using 1:2 with lines, 'out.rk4' using 1:2
Notice that the anhbust integration is tracking the rk4 solution quite well, but with much bigger jumps.

For a different perspective, try a phase space plot of y vs dy/dt with 

   plot 'out.bust' using 2:3 with lines, 'out.rk4' using 2:3