ODE: My Ordinary Differential Equation Program

Andrew D. Lewis

18 March 1998

1 Introduction

This is documentation I wrote for my ordinary differential equation solver which I call ODE. The documentation was written mainly so I would remember how to use the program after long periods of enforced absence.

The program is set up to allow easy parameter changes, and convenient control of the output. The program also produces a log file which should allow the user to reproduce any data at a later date.

The numerical method is simply a fourth order Runge-Kutta scheme. I may change this later to be something cooler.

2 Program Compilation^{1}

To get a version of ODE up and running, simply edit the variables BINDIR and MANDIR in Makefile to reflect your system settings. Also, DOCDIR and DEMODIR are used to build a correct man page. Then type make to build the binary, make install to put it in the specified location. make man will build and install the man page. make clean will get rid of object files, etc.

3 Basic Usage

To use ODE, simply type

ode filename[.dat]

If the argument is left off, ODE will read input from standard input and write output to standard output. If the argument is present, ODE will check for the existence of a file filename.odeinp from which to read the input data. If this file is not present, ODE will next check for a file odeinp.odeinp. If this file is also not present, you will have to enter input data by hand. The input data consists of the following things:

1This set up procedure has been performed with SunOS 5.4 and is not guaranteed to work on any other operating system.

1. The parameters for the problem. In Section 5 we will tell you how to set up ODE to give you the parameters you need.

2. Next you enter the initial conditions.

3. Next enter the beginning time for the simulation.

4. Next enter the time for the simulation to end. You may run simulations backwards.

5. Next enter the time step.

6. Next you specify the number of degrees of freedom you wish to print to the output file.

7. Finally you specify which degrees of freedom you wish to print out.

Now ODE will take this data and integrate the differential equation. The output data will be place in the file filename.dat. A log file, filename.odelog will be created which contains the time of the simulation, as well as all the input data. If the output is written to standard output, no log file is created.

4 Specifying Which Differential Equation will be Integrated

Now you know how to make ODE integrate an ordinary differential equation. But which differential equation will it integrate? Well, the one you define with the subroutine derivs. In the src directory you will see a file derivs.c which contains the definition of the ordinary differential equation you will integrate. Every time you change this file, you must recompile ODE to reflect these changes. The number of degrees of freedom is specified in ode.h by the variable NDOF. Any changes here will also necessitate a recompilation of ODE.

The format of derivs.c is fairly simple. The vector y is the state vector. The vector ydot contains the right hand side vector for the differential equation. The vector parms contains all parameters. Now type…

Again, let me restate that you must change NDOF in ode.h to specify the correct number of degrees of freedom.

5 Specifying the Parameters in the Differential Equation

You will need to specify the values in the parms vector in derivs.c. This is done in the subroutine parmenter which is in the file parmenter.c. This subroutine should be set up so that it echoes to the screen which parameter is to be read, and then allows the data to be read from either the screen or from the input data file. If you wish, you may want to echo the value of the parameter to the log file.

The following segment of code illustrates how one may read in parms[0].

    printf("\nEnter the value of mass in kg.\n");
    fscanf(input_fp,"%lf",&m);
    if (!UsingStdIn) {
      fprintf(log_fp,"\nThe value of mass is: %lf kg\n",m);
    }

I hope this is self-explanatory. I should also mention that you may wish to use parmenter.c to do some computations involving parameters that you would otherwise have to do by hand. This might allow you to, for example, enter a rotational speed in rpm, but convert it to rad/sec for use in the program. I have found this very helpful.

6 Initial Condition Hacks

I often am in a situation where I need some of my initial conditions to be determined by the values of others. To facilitate this, I use the subroutine ichacks which is in the file ichacks.c. This subroutine simply specifies some values of a vector y according to rules you make. For example, if you always wanted y[2] to be half the value of y[0] at the initial time, you would write

    y[2] = y[0] / 2.0;

in ichacks.c. ichacks is always called after the input data is read in. Thus, in the above case, the initial value of y[2] which you specify when you enter the input data would be irrelevant since it would be overwritten by ichacks.

7 Demo

There is a demo set up in the demo directory. It will integrate a mass-spring system, and should give you some idea of how the subroutines derivs, parmenter, and ichacks work.

Please send any comments and errors to: andrew@mast.queensu.ca.