PHYCS 6730 Lab Exercise: Ordinary Differential Equations: Step Size Error

In these exercises we will be comparing the error as a function of step size for the modified Euler and 4th order Runge Kutta methods for integrating an ordinary differential equation. The code is in ~p6730/exercises/ode and in the course library directory ~p6730/examples/nmr. We will be comparing results from checkeul2 and checkrk4.

Exercise 1 The Code

Please look at the codes checkeuler2.cc and checkrk4.cc. These codes do a single modified Euler step and a single 4th order Runge Kutta step, respectively.

Notice that both codes are solving a pair of first order equations. Throughout the code the dependent variables are kept in a two-dimensional array y and their derivatives, in dydt. Think of y[0] as an angle theta and y[1] as an angular velocity omega. When you run the code, it prompts for three values, theta0, omega0, and dt. These are the values of theta = y[0], omega = y[1] at time t = 0 and the integration time step, which is called h in the lecture notes .

After looking over the code, please answer these questions:

  1. What are the two first order differential equations being integrated? How are they reexpressed as a single second order differential equation? Hint: The derivatives are evaulated in the subprogram derivs.
  2. Why is it necessary to call derivs before calling NR::rk4 and euler2? Hint: Check the conventions (and instructions) for NR::rk4 and euler2. The code for euler2 is in ~/examples/nmr/euler2.cc.
  3. Check that the exact solution to the differential equation for initial conditions y[0] = 1 and y[1] = dy[0]/dt = 1 is y[0] = cos(t) + sin(t). Do this by direct substitution.
  4. In terms of the parameters that the program asks the user to supply, what are the starting and ending times for this integration? Hint: What is the initial time, the step size, and number of steps?

Exercise 2 Experimenting with local step size error

Run each code with initial theta = 1 and omega = d theta/dt = 1 and dt chosen from the list below. Complete the table 
    dt     y_exact(dt)  y_Euler2(dt)    diff       y_rk4(dt)     diff
    .1     1.094837582  1.09500000  -.000162418    1.09483750   .82e-7
    .2     
    .4
Verify that the error in the modified Euler result (in a single step) varies as dt^3 and in the 4th order Runge-Kutta as dt^5. To show this you should divide the difference values in the table by the appropriate power of dt and see if the ratios look like they are approaching a constant as dt goes to zero (as opposed to blowing up if the power was really lower or vanishing if the power was really higher.)