# 780.20: 2082 Session 6

Handouts: diffeq_test.cpp and diffeq_routines.cpp, excerpt from Chap. 9 of "Computational Physics".

In the chapter 9 handout is an introduction to algorithms for integrating differential equations. In this session, we'll go over the basic ideas by example, using the routines in session06.tarz, in preparation for looking at "Anharmonic Oscillations" and "Differential Chaos in Phase Space".

• Run a code to integrate a simple first-order differential equation using Euler's and 4th-order Runge-Kutta algorithms, then modify it so you can analyze the errors.
• Add code for a 2nd-order Runge-Kutta algorithm.
• Adapt the code to treat the 2nd order F=ma problem in chapter 9 of the text.

Please work in pairs (more or less). Dick and Daniel will bounce around 2082 and answer questions.

## Integrating a First-Order Differential Equation

The code diffeq_test.cpp calls the differential-equation-solving routines in diffeq_routines.cpp ("euler" and "runge4") to integrate a series of coupled differential equations (but we'll start with a single equation). Functions for the right-hand-side of a first order differential equation (dy/dt = rhs[y,t]) and the exact y(t) [called "exact_answer"] for a specified initial condition are also defined in this file. There is also a header file diffeq_routines.h.

1. Look through the Session 6 Background Notes for a quick overview of differential equation solving.
3. Use make_diffeq_test to compile and link diffeq_test. Run the program to generate diffeq_test.out and look at it in an editor. The gnuplot plotfile diffeq_test.plt generates comparison plots of the integrated function from the output in diffeq_test.out. Load this plotfile in gnuplot:
and examine the result. What can you conclude at this point?

4. Look at the printout for diffeq_test.cpp and diffeq_routines.cpp and compare to the handout with the excerpt from Chapter 9 of Landau and Paez to figure out what is going on. The codes follow the notation in the excerpt, so t is the independent variable and y[] is an array of solutions (see section 9.5 in the handout). At present there is only one equation (first-order), so only y[0] is used. The routines euler and runge4 are based on Eqs. (9.34) and (9.46), respectively. What is the differential equation being integrated?

5. Modify the diffeq_test.plt file to plot the relative error at each value of t. (Modify the plot file and NOT the program; see the gnuplot handout on plot files for an example of how to do this.) As usual in studying errors, a log-log scale will be useful. The first point at t=0 may get in the way. Use "xrange = [?:?]" in gnuplot (where you fill in the ?'s) to avoid this problem. What can you say qualitatively about the errors?

6. Now generate and plot results for a second value of h (your plot should have both values of h, so think about how to best do this). You'll want it use something like 1/10 the value, so it's easy to check the effect (e.g., if the difference goes like hn, then you'll see 10n, which is easy to see on a log-log plot). When the local error (for each step h) adds coherently, then the "accumulated" or "global" error for a given algorithm at tf scales like Nf*(local error) = (tf/h)*(local error). You can verify for Euler's algorithm, for which the local error to be h2 (see notes and excerpt), that the global error is, in fact, h. What is the local error for 4th-order Runge-Kutta according to the graph?
7. Next, check how the accumulated error at one value of t scales with h for the two algorithms. Take t=1, for example. You'll need to modify the code to output the results for y(t=1) for a range of h's (think logarithmic!). Some things to be careful of:
• Print out enough digits. For small h's, 9 digits is not enough (try 16).
• The most common problem is printing out y_rk4[0] and exact_answer(t,params_ptr) at two different points. (Note: it is not important that the t used is exactly the same for every h, but it must be the same for the exact and rk4 result for each h.) Look at your output file!
There should be different regions of the error plot, as we've seen before. Interpret them.

Given your results, how would you choose a step size for 4th-order Runge-Kutta? (Hint: How do you explain the behavior of the error for small h?)

## 2nd Order Runge-Kutta Algorithm

This exercise is intended to verify that you understand the meaning of the Runge-Kutta algorithms and how they are implemented in C++.

1. Add a third diffeq routine to the code, which implements the 2nd-order Runge-Kutta algorithms described in section 9.9 [and summarized in Eqs. (9.44) and (9.45)].
2. Check the scaling of the error with h.
3. Try a different differential equation, such as dy/dt = 2*cos(2pi*t).

## Integrating a 2nd-Order Differential Equation

Second-order differential equations are treated by writing them as two coupled 1st-order equations, as described in Landau and Paez, Chapter 9. We'll try this out for a simple example, which we'll generalize later to look at chaotic behavior.

1. Consider the differential equation for a simple, undriven harmonic oscillator [e.g., a mass m on a spring with constant k: d2x/dt2 + (k/m)x = 0]. What is the general solution in terms of the initial position and velocity?

2. Rewrite the differential equation as two coupled 1st-order equations [for x and v=dx/dt]. Use units in which the oscillator mass and spring constant are equal to one. Take the initial position to be 1.0 and the initial velocity to be zero.
3. Modify a copy of diffeq_test.cpp to use the runge4 subroutine to calculate this oscillator for times t=0 to t=10. You'll need to change N, modify rhs to consider both i==0 and i==1, put the exact answer in exact_answer, and change the limits of t appropriately.
4. Plot the result and the exact result for comparison. What do you conclude?