//////////////////////////////////////////////////////////////////////////////// // 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; }