A large number of man-made systems and natural phenomena can be conveniently described by using differential equations. These equations generally involve derivatives and integrals of the dependent variables with respect to one independent variable. The graph of the solution of the differential equation describing a system is often useful for analysis of the process.
In this chapter we will develop numerical methods for solving these differential equations and displaying the result in graphical form. We start with developing a method to solve first order differential equations numerically using the fourth order Runge-Kutta method. The fourth order Runge-Kutta method is one of the standard (perhaps the standard) algorithm for numerical solution of differential equations. It is generally considered to provide an excellent balance of power, precision and simplicity to program. This method can be extended to solve second and higher order differential equations.
The Runge-Kutta Method
In this section we consider the first-order differential equation that relates a pair of variables. Such equations may be written in the form:
y'(t) = g(t, y)
where the ' indicates differentiation with respect to t.
The fourth order Runge-Kutta is based on a fourth order approximation which requires four evaluations of the function g(t, y) to determine the next value of the dependent variable y(t). The four terms, gn, below are defined at some value ti of the independent variable t as follows:
g1 = g(ti, yi)
g2 = g(ti + dt/2, yi + g1·dt/2)
g3 = g(ti + dt/2, yi + g2·dt/2)
g4 = g(ti + dt, yi + g3·dt)
where:
dt = ti+1 - ti
The expressions above must be calculated in the order given since g2 depends on the value of g1, etc. After the terms above have been calculated, the next value of y(t) is found using the relation:
yi+1 = yi + dt·(g1 + 2·g2 + 2·g3 + g4)/6
This is the fourth-order Runge-Kutta method; it is easily programmed in C.
Coding the function: g(t, y)
This function may take any of several forms; continuous, piecewise-linear, or discontinuous. In order to accommodate as many forms as possible we will define a general C function that returns a double and takes two double arguments, the first for the independent variable t and the second for the dependent variable y. The code within the function will implement the right hand side of the differential equation. For example, to implement the differential equation:
y' = -y
The function G could be encoded as:
double G( double t, double y )
{
return ( - y );
}
Coding the Runge-Kutta algorithm
An algorithm that computes the numerical solution to a first order differential equation must be given the following information:
It must then return a value of the dependent variable y(t) at the end of a computation interval.
A pseudocode representation of the process is shown here:
begin
dt := computation interval / segments
t := initial value
y := initial value
i := 0
loop
exit when i = segments
g1 := g(t, y)
g2 := g(t + dt/2, y + dt*g1/2)
g3 := g(t + dt/2, y + dt*g2/2)
g4 := g(t + dt, y + dt*g3)
y := y + dt*(g1 + 2*g2 + 2*g3 + g4)/6
t := t + dt
i := i + 1
end loop
end
A function in C that will perform the Runge-Kutta algorithm is as follows:
double rk( double t_beg, // interval start
double y_beg, // derivative at start
double delta, // interval length
int n, // number of segments > 0
double (g)(double, double) ) // integrand function
{
double g1, g2, g3, g4;
double t = t_beg; // initialize time
double y = y_beg;
double dt = delta / n; // initialize segment period
int i; // loop counter
for( i = 0; i < n; i++ )
{
g1 = g( t, y );
g2 = g( t + dt / 2.0, y + dt * g1 / 2.0 );
g3 = g( t + dt / 2.0, y + dt * g2 / 2.0 );
g4 = g( t + dt, y + dt * g3 );
y += ( g1 + 2.0 * g2 + 2.0 * g3 + g4 ) * dt / 6.0;
t += dt;
}
return y;
}
Using the Runge-Kutta algorithm
In order to draw a curve representing a solution to a differential equation we need to include the functions above in a program loop that can calculate points through which a curve can be drawn. In general the computation should do the following:
t := t0
y := y0
points := number of points for defining a curve
delta := step size in units of t
a0 := plot_point(t, y)
i := 0
loop
exit when i = points
y := rk(t, y, delta, segments, g(t, y))
t := t + delta;
ai := plot_point(t, y)
end loop
A C function for solving the differential equation:
y' = -y, y(0) = 1
can be coded as follows:
int calculate( POINT *a0, int points )
{
int i;
double y = 1.0;
double t = 0.0;
double delta = 0.05;
a0[0] = np.plot_point( t, y );
for( i = 1; i < points; i++ )
{
// next iteration
y = rk( t, y, delta, 2, G );
t = t + delta;
// append results to curves
a0[i] = np.plot_point( t, y );
}
return i; // items in POINT array
}
in which the function G(t, y) is as shown above.
Higher order differential equations
Many processes must be represented by more complex differential equations. For example, a third order differential equation:
y''' + 2·y'' + 2·y' + y = u
where y(t) is the dependant variable to be solved numerically, u(t) is the forcing function, and t is the independant variable.
In general, an nth order differential equation can be decomposed into a set of n first order differential equations. For the third order differential equation above, define:
x1 = y(t)
x2 = y'(t)
x3 = y''(t)
Then the third order differential equation can be written as three first order equations:
x'1 = x2 x'2 = x3 x'3 = - x1 - 2·x2 - 2·x3 + u
These first order state equations can be represented in vector - matrix form. Letting:
x1
x = x2
x3
Then, the state equations are:
x' = A x + b u
Where,
0 1 0 0
A = 0 0 1 b = 0
-1 -2 -2 1
The vector - matrix form can also be applied to non-linear state equations; in general:
x'(t) = f( x(t), u(t), t )
Example
This example uses the Runge-Kutta numerical differentiation method to solve the matrix differential equation described above.
/*----------------------------------------------------------------------------*\
** Put your mathematical functions to plot here. **
\*----------------------------------------------------------------------------*/
int calculate( POINT *, POINT *, POINT *, int );
vector mat_rk( double, vector &, double, int, vector (g)( double, vector & ) );
vector G( double t, vector & );
/*----------------------------------------------------------------------------*
** "run" performs calculations and plots results.
**
** This is a shell in which the function calls for a specific calculation
** are entered.
**
** returns: nothing
*/
void run()
{
POINT curve0[110], curve1[110], curve2[110];
int i;
// first, define a plot region
// this must be done before calculating points to be plotted
np.plot_region( -0.1, 10.0, 12, -3 );
// do your calculations
i = calculate( curve0, curve1, curve2, 110 );
// plot the results
np.draw_axes();
np.red_curve( curve0, i );
np.green_curve( curve1, i );
np.blue_curve( curve2, i );
}
/*----------------------------------------------------------------------------*
** "calculate" is a shell for the calculation algorithms.
**
** returns: number of items in POINT array
*/
int calculate( POINT *a0, // curve 0
POINT *a1, // curve 1
POINT *a2, // curve 2
int points ) // count
{
int i;
vector y;
double t = 0.0;
double delta = 0.1;
a0[0] = np.plot_point( 0, 0 );
a1[0] = np.plot_point( 0, 0 );
a2[0] = np.plot_point( 0, 0 );
for( i = 1; i < points; i++ )
{
// next iteration
y = mat_rk( t, y, delta, 2, G );
t += delta;
// append results to curves
a0[i] = np.plot_point( t, y[0] );
a1[i] = np.plot_point( t, y[1] );
a2[i] = np.plot_point( t, y[2] );
}
return i;
}
/*----------------------------------------------------------------------------*
** "mat_rk" performs numerical differentiation using the 4-th order
** Runge Kutta algorithm to solve matrix differential equations of the form:
**
** y' = g( t, y )
**
** returns: nothing
*/
vector mat_rk( double t_beg, // interval start
vector &y_beg, // y at start
double delta, // interval duration
int n, // number of steps > 0
vector (g)( double, vector & ) ) // integrand function
{
vector y = y_beg; // initialize vector
double t = t_beg; // initialize time
double dt = delta / n; // initialize step duration
vector g1, g2, g3, g4; // variables for computation
double dt2 = dt / 2.0;
double dt6 = dt / 6.0;
int i; // loop counter
for( i = 0; i < n; i++ ) // do n steps of R-K algorithm
{
g1 = g( t, y );
g2 = g( t + dt2, y + g1 * dt2 );
g3 = g( t + dt2, y + g2 * dt2 );
g4 = g( t + dt, y + g3 * dt );
y = y + ( g1 + g2 * 2.0 + g3 * 2.0 + g4 ) * dt6;
t = t + dt;
}
return y;
}
/*----------------------------------------------------------------------------*
** "G" defines a matrix differential equation of the form:
**
** y' = G(t) y + u(t)
**
** returns: the left hand side of the above
*/
vector G( double t, vector &y )
{
vector g, u;
// the input vector u(t)
u[0] = 0;
u[1] = 0;
u[2] = 10;
// the equation
g[0] = 0.0 * y[0] + 1.0 * y[1] + 0.0 * y[2] + u[0];
g[1] = 0.0 * y[0] + 0.0 * y[1] + 1.0 * y[2] + u[1];
g[2] = - 1.0 * y[0] - 2.0 * y[1] - 2.0 * y[2] + u[2];
return g;
}