Solving ODEs in C

I set my students an exercise: to solve the nonlinear differential equation

\ddot{x} + \sin(x) = A\sin(\omega t).

What follows is an outline of a model solution (but by no means complete, or optimal). I take this opportunity to introduce several other features of C that I haven’t had the chance to talk about in the classes.

This equation is, of course, the (nondimensionalised) equation of motion of a forced pendulum, where x represents the angle that the pendulum makes with the downward vertical. If A = 0 or \omega = 0, then the pendulum is unforced, and executes regular oscillations. The period of these oscillations is approximately equal to 2\pi for small-amplitude oscillations (where the motion is approximately simple harmonic motion), but increases as the amplitude increases. If forcing is turned on, then a second period—that of the forcing—is introduced. The pendulum may still oscillate with its ‘intrinsic’ period, or if the forcing is strong enough then the pendulum may start oscillate with the period of the forcing. If the period of the forcing is close to the intrinsic period of the pendulum, then resonance may occur.

We first write the equation as a system of two equations,

\dot{x} = y \text{ and } \dot{y} = - \sin(x) + \sin(\omega t)

which is first-order. We then solve this system numerically using the Runge-Kutta method. The program diffsolve is meant to do this using a Runge-Kutta solver. Its full source code can be downloaded here. The code, as it is written, is quite general; it can solve any system of two ODEs of the form

\dot{x}=f(t,x,y) \text{ and } \dot{y}=g(t,x,y).

Annotation of the code

After the initial comments and the #includes, we first come across a typedef struct:

typedef struct {
    double amp;
    double omega;
} problem_t;

I declare a type, problem_t, which will store the parameters of the problem, amp and omega. I usually wrap all parameters into a struct like this, because it lets me pass a single variable, of type problem_t, to a function, instead of passing each of the individual parameters. It also means that, if I want to change my program to allow more parameters, I can simply change the declaration of problem_t, rather than changing the definition of every function that needs these parameters. This pays off when we work with problems with ten or twenty parameters.

The first function defined makes use of this:

double forcing(double t, problem_t params) 
{
    return params.amp * sin( params.omega * t );
}

This is used by the function rhs_y, which calculates \dot{y} . (You could do away with the function forcing and have the above inside rhs_y, if you prefer.)

The function rhs_x, which gives the derivative \dot{x} , appears to require arguments that it does not use:

double rhs_x(double t, double x, double y, problem_t params)
{
    return y;
}

Although the arguments t, x and params are not actually used, I have chosen to include them anyway, because in the future I might want to move from solving \dot{x} = y to something more general, of the form \dot{x} = f(t,x,y) .

Now we come to the main function. The first thing we notice is that the function prototype says

int main(int argc, char* argv[]) 

instead of simply

int main()

We want to be able to run our program without having to recompile the program each time we want to change our parameters such as initial conditions. We will type something like

$ ./diffsolve 1 2 100 0.005 1 0

into the command line, in order to specify the different parameters to diffsolve. These are command line arguments, and they are treated by C as arguments to main. The argument argc is the argument count, the number of arguments; and argv is the argument vector, an array of strings. Note that the program’s name, ./diffsolve, is treated as the 0th argument, and so argc is always at least equal to 1.

The arguments are actually read in the following block:

if (argc > 6) 
{
    params.amp   = atof(argv[1]);
    params.omega = atof(argv[2]);
    tmax  = atof(argv[3]);
    dt    = atof(argv[4]);
    x     = atof(argv[5]);
    y     = atof(argv[6]);
}
else
{
    fprintf(stderr, "Usage: %s amp omega tmax dt x0 y0\n", argv[0]);
    exit(-1);
}

If the user has specified at least seven arguments (where the 0th argument is the name of the program), then these arguments are taken. The atof function converts strings to floats (although our variables are declared as doubles; this is fine, because floats can be implicitly cast into doubles).

If the user has not specified enough arguments, then the program gives an error which tells the user how to use it, and exits with a return code of -1, which (by convention) indicates an error. (The exit function is provided in stdlib.h.) The error message is printed using fprintf, which stands for ‘printf to a file’, or ‘print formatted text to a file’. Here the ‘file’ is stderr, the standard error stream. In contrast, the function printf writes to the standard output stream, represented by stdout. By default, both stdout and stderr write to the screen, but the two streams may be separated from each other. It is a good idea to use stderr for errors, warnings or diagnostic messages, and stdout for actual output that you want to process.

The rest of the program is fairly straightforward: the values of t, x and y are updated in a while loop, and at each time step their values are printed (to standard output).

Usage (on Unix-like systems)

After compiling, the command

$ ./diffsolve 1 2 100 0.005 1 0

will print plenty of lines, giving the solution of the system at each timestep.

You probably want to save this information to another file, which you can do by writing

$ ./diffsolve 1 2 100 0.005 1 0 > out.txt

This will write the results to the file out.txt, which you can then process somewhere else. For example, you can use gnuplot to plot the results:

$ gnuplot
gnuplot> set terminal png
gnuplot> set output "out.png"
gnuplot> plot "out.txt" using 1:2 with lines

This produces the image out.png:
out.png

Leave a Reply

Your email address will not be published. Required fields are marked *