Modular Programming

We programmed in previous sections the Euler and Mid-Point methods for a scalar (1-dimensional) and simple mechanics examples from within a simple main code body. However, good programming practice demands that the details be relegated to functions. This allows the overall logic of the solution to be clearly exhibited in the main program.

The code must be able to solve simultaneously for many unknown functions (vector of functions) in the case of coupled ODEs.

Furthermore, in the case of solving ODEs, much is to be said for a clean separation of the methods of solutions from the details of the model. A properly written code should allow adapting the solution of a model to another with minimum effort. Also the same code should allow the use of another method of solution to be trivial, i.e., a switch from a Euler solution to a Mid-Point solution should be easy.

This argues for interchangeable C functions which implement a time step within the context of whatever method, e.g., a Euler step or a Mid-Point step. This function should input the current time, the current time step, and the current function values. A prototype could be

void mid_point_step ( double t, double dt, double f_t[] )

A general system of ODEs implies one or many unknown functions. The implementation in C requires arrays of specific dimension corresponding to the dimension of the space in which the solution leaves.

The time step function requires the availability of a C function that calculates the Right Hand Sides of the differential equations, themselves contained in an array. This function should accept as input the current time, the current function values (in an array), and output the RHS values (in an array). A prototype could be

void rhs_function( double t, double f[], double rhsf[] )

The name of this function can be set ahead of time or be provided to the time step function as a function pointer.

A clear advantage of this modular approach is that the time step function can be concisely programmed in a manner that closely follows the mathematical statement of the algorithm. For instance, the Mid-Point time step function could be written as mid_point_step.c The Euler time step function could look like euler_step.c Note that all the functions must be marched simultaneously in time.

The code HO_euler.c implementing the solution of a mass on a spring serves to illustrate the use of the time step routine. The model dimension and parameters are defined as global variables. This allows the time step routine to set temporary arrays of appropriate dimensions for its inner workings. The model parameters are only needed by the RHS function (and the conserved total energy).

The Euler and Mid-Point integrator step functions using function pointers are euler_step_pointer.c and mid_point_step_pointer.c

HO_euler_pointer.c implements function pointers.