////////////////////////////////////////////////////////////////////////////////
// File: graggs_method.c //
// Routines: //
// Graggs_Method //
// Graggs_Method_Richardson //
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// //
// Description: //
// Gragg's method, also called the modified midpoint method, for approxi- //
// mating the solution of the differential equation y'(x) = f(x,y) with //
// initial condition y(x0)= c is a second order method based upon the //
// following theorem due to Gragg. //
// //
// Thm: Let f(x,y) have continuous partial derivatives up to order 2N+2, //
// y(x) be the exact solution to the differential equation y' = f(x,y) //
// with y(x0) = y0, let x[i] = x0 + i*h, and let Y(x,h) be defined //
// inductively by: Y(x0,h) = y0, Y(x[1],h) = y0 + h*f(x0,y0), and //
// Y(x[i],h) = Y(x[i-2],h) + 2h*f(x[i],Y(x[i],h)) for i >= 2. Then //
// Y(x,h) = y(x) + Sum(h^2k [uk(x) + (-1)^(x-x0)/n vk(x)] + h^(2N+2)e(x,h)//
// where the Sum extends from k = 1 to N, for all x in the domain of //
// definition, all h = (x - x0)/n, n = 1,.... The function e(x,h) for //
// fixed x is a bounded function of h. //
// //
// The oscillating term (-1)^(x-x0)/n vk(x) can be removed by noticing //
// that Yest(x,h) = 1/2[Y(x,h) + Y(x-h,h) + h*f(x,Y(x))] //
// has an expansion using the theorem above //
// Yest(x,h) = y(x) + h^2 [u1(x) + 0.25*y''(x)] + ... //
// in which the oscillating term no longer appears in the leading term. //
// //
// Richardson extrapolation can be used to increase the order and //
// accuracy. //
// //
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// double Graggs_Method( double (*f)(double, double), double y0, double x0, //
// double x, int number_of_steps ); //
// //
// Description: //
// This function solves the differential equation y'=f(x,y) with the //
// initial condition y=y0 at x = x0. The return value is the value of //
// the solution y evaluated at x. //
// //
// Arguments: //
// double *f //
// Pointer to the function which returns the slope at (x,y) of the //
// integral curve of the differential equation y' = f(x,y) which //
// which passes through the point (x0,y0). //
// double y0 //
// The initial value of y at x = x0. //
// double x0 //
// The initial value of x. //
// double x Final value of x. //
// int number_of_steps //
// A positive even number, the initial step size is //
// is h = (x - x0) / number_of_steps. //
// //
// Return Values: //
// The solution y(x) of the initial value problem y' = f(x,y) with //
// y(x0) = y0. //
// //
////////////////////////////////////////////////////////////////////////////////
// //
double Graggs_Method( double (*f)(double, double), double y0, double x0,
double x, int number_of_steps ) {
double h = (x - x0) / (double) number_of_steps;
double h2 = h + h;
double y1 = y0 + h * (*f)(x0,y0);
double y2;
while ( --number_of_steps ) {
x0 += h;
y2 = y0 + h2 * (*f)(x0,y1);
y0 = y1;
y1 = y2;
}
return 0.5 * ( y0 + y1 + h * (*f)(x,y1) );
}
////////////////////////////////////////////////////////////////////////////////
// double Graggs_Method_Richardson( double (*f)(double, double), double y0, //
// double x0, double x, int number_of_steps, int richardson_columns ) //
// //
// Description: //
// This routine uses Gragg's method together with Richardson extrapolation//
// to solve the differential equation y'=f(x,y) with the initial condition//
// y = y0 at x = x0. The return value is the value of the solution y //
// evaluated at x. //
// //
// Arguments: //
// double *f //
// Pointer to the function which returns the slope at (x,y) of the //
// integral curve of the differential equation y' = f(x,y) which //
// which passes through the point (x0,y0). //
// double y0 //
// The initial value of y at x = x0. //
// double x0 //
// The initial value of x. //
// double x Final value of x. //
// int number_of_steps //
// A positive even number, the initial step size is //
// is h = (x - x0) / number_of_steps. //
// int richardson_columns //
// The maximum number of columns to use in the Richardson //
// extrapolation to the limit. //
// //
// Return Values: //
// The solution y(x) of the initial value problem y' = f(x,y) with //
// y(x0) = y0. //
// //
////////////////////////////////////////////////////////////////////////////////
// //
static const double richardson[] = {
1.0 / 3.0, 1.0 / 15.0, 1.0 / 63.0, 1.0 / 255.0, 1.0 / 1023.0, 1.0 / 4095.0
};
#define MAX_COLUMNS 1+sizeof(richardson)/sizeof(richardson[0])
#define max(x,y) ( (x) < (y) ? (y) : (x) )
#define min(x,y) ( (x) < (y) ? (x) : (y) )
// //
////////////////////////////////////////////////////////////////////////////////
// //
double Graggs_Method_Richardson( double (*f)(double, double), double y0,
double x0, double x, int number_of_steps, int richardson_columns ) {
double dt[MAX_COLUMNS]; // dt[i] is the last element in column i.
double integral, delta;
int j,k;
richardson_columns = max(1, min(MAX_COLUMNS, richardson_columns));
for (j = 0; j < richardson_columns; j++) {
integral = Graggs_Method( f, y0, x0, x, number_of_steps );
for ( k = 0; k < j; k++) {
delta = integral - dt[k];
dt[k] = integral;
integral += richardson[k] * delta;
}
dt[j] = integral;
number_of_steps += number_of_steps;
}
return integral;
}