The Simplest ODE System of Equations

Consider the differential equation

y'' + y = 0

where y' means the first derivative of y, y'' is the second derivative of y.

We decompose this equation into two first-order equations with the following method:

y1 = y
y2 = y'

Now there will be two coupled differential equations, let's see,

y1' = y2
y2' = -y1

This is the simplest ODE system of equations because we already know the answer, so it is an easy matter to check the numerical solution. The solution is a sinusoid of type that depends on the initial conditions.

For this problem we choose an equally simple, reliable second-order numerical integrator called "The trapezoid rule," or "Heun's method," or "the implicit Euler method," or "the Backward Difference method." For example, for a single ODE this would be

y(n+1) = y(n) + h * f(x(n), y(n))

the predictor, and

y(n+1) = y(n) + (h/2) * (f(x(n), y(n)) + f(x(n+1), y(n+1)))

the corrector.

We wish to integrate this with the new ODE systems integrator in HMD. It can easily be compared to an integrator hard-coded in C or Fortran.

The fundamental part of the ODE is the derivative function, y' = f(x,y). Or for the system of ODE is is the two functions

f1(x, y1, y2)
f2(x, y1, y2)

In HMD these are callback functions. First define the functions in an ascii file. From the HMD interpreter you will import these with the source command. Or just put all the script in an ascii file and run it from the command line line,

hmd -f ode1.hmd

//
//  The derivative function for y1
//
function fxy1 { double x, double y1, double y2 } {
    double ret;

    ret = y2;

    return ret;
}

//
//  The derivative function for y2
//
function fxy2 { double x, double y1, double y2 } {
    double ret;

    ret = -1 * k * y1;

    return ret;
}

Each time the ODE stepper is called in a loop the stepper will "call back" these two functions. The HMD platform will supply as arguments the current position, x, and the state variables y1, y2.

Now, to the script command that set up the ODE machinery. There is a prepackaged Heun integrator supplied with HMD. You don't need to worry about entering the weight tables and other parameters. Both equations will be set up as Heun models with the following commands. Given that we define the initial conditions and step parameters with

hx = 0.1;
x0 = 0.0;
y0 = 1.0;       // init state 1
dy0 = 0.0;      // init state 2
N = 63;         // number of steps
k = 1.0;        //  constant in y'' + k*y = 0

The we begin setting up the models with the script functions,

h_odemodel1 = ode_heun_open(hx, x0, y0, 0, fxy1, 3, $1, $2, $3);

h_odemodel2 = ode_heun_open(hx, x0, dy0, 0, fxy2, 3, $1, $2, $3);

Notice that the ode_heun_open functions take the names of the callback functions as parameters: fxy1, and fxy2. The argument following the callback name is an integer specifying how many placeholder symbols follow. The placeholder symbols represent the independent variable, x, and the two state variables y1 and y2. For the moment, these three placeholders are merely telling HMD that there will be three parameters to the callback functions.

Now we need to tell HMD that these two models (ODE equations) will be part of a system of equations. The function for initializing the system is init_ode_system, which will be called as

h_odesystem = init_ode_system(h_odemodel1, h_odemodel2);

The function takes as parameters the handles for the individual ODE models that were created above. It returns a new handle that identifies the system of equations.

Because we are creating a system of coupled equations, not just some single ODE equations, an extra specification of the callback arguments is required. By default, for a single ODE, the symbolic placeholders mean $1 for x, $2 for y, and higher numbers mean higher derivatives of y. But for a system of equations we expect to pass multiple state variables to any of the integrator callbacks. This is done with the set_function_args_ode_system call.

set_function_args_ode_system(h_odesystem, h_odemodel1, 1, 1, '$, $1_0, $2_0');

set_function_args_ode_system(h_odesystem, h_odemodel2, 1, 1, '$, $1_0, $2_0');

In the above function calls the callback specifiers are given within strings defined by single quotes. Here $ means x, $1_0 means the zeroth derivative of the first state variable, and $2_0 means the zeroth derivative of the second state variable. Compare this to the function definitions of fxy1 and fxy2 shown above.

Now before starting the integration it would be nice to have the output data saved to arrays. So we will "attach" data tables to each of the ODE models. If a model has been given a data table it will automatically be written with data during the integration. On the other hand, if you never attach a data table the model can be integrated, but no data would be saved. That is probably not what you want. So attaching data tables is performed as follows,

add_datatable_ode(h_odemodel1, N);
add_datatable_ode(h_odemodel2, N);

Notice that we make the data table the same size as the number of steps with the parameter N.

Before starting the integration there are a couple of last initializations to be performed. These are

update_states_ode_system(h_odesystem);

and

compute_derivatives_ode_system(h_odesystem);

The integration of the ODE or system of ODEs centers on the functions step_ode and step_ode_system. For a system of equations we use the latter.

step_ode_system(h_odesystem, 'repeat', N);

Now to get the output data. We request the data from the ODE system object by using the handle to the object.

outstates = get_ode_system_data(h_odesystem);

which you give you a struct in the HMD workspace. This struct has members rows, the number of rows in the data matrix (in this case 3), columns, the number of columns in the matrix which is also the number of integration steps, and data being the actual array of data. Our array, outstates.data, is a binary array in the same sense as a C language block of data. You can view this data in the HMD workspace or save it to a file using the dump_block function. Let's save the data to a file called "odesystem1.dat," which will be an ASCII file with three columns: the x data, the y1 data, and the y2 data.

dump_block(outstates.data, 1, 'odesystem1.dat');

Now clean up the HMD workspace,

ode_close(h_odemodel1);
ode_close(h_odemodel2);
close_ode_system(h_odesystem);

The reason why this is the simplest ODE system of equations is because we can readily check the truncation error of the answer. The extrema of the solution should be at -1 and +1.