//////////////////////////////////////////////////////////////////////////////// // File: bulirsch_stoer.c // // Routines: // // Gragg_Bulirsch_Stoer.c // //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// // // // Description: // // The Gragg-Bulirsch-Stoer method for approximating the solution of the // // differential equation y'(x) = f(x,y) with initial condition y = c when // // x = x0 is a variable step method which uses Graggs modified midpoint // // method for a given step size h to evaluate y(x,h), a rational function // // approximation to the limit a h -> 0, and an algorithm to estimate the // // error and adjust the step size so that the error is within a // // preassigned tolerance. // // // //////////////////////////////////////////////////////////////////////////////// #includestatic int number_of_steps[] = { 2,4,6,8,12,16,24,32,48,64,96,128 }; #define ATTEMPTS sizeof(number_of_steps)/sizeof(number_of_steps[0]) #define LENGTH_SCALE 1+(ATTEMPTS>>1) static double Graggs_Method( double (*f)(double, double), double y0, double x0, double x, int number_of_steps ); static int Rational_Function_Zero_Extrapolate( double *zero, double x[], double f, double aux_array[], int n ); //////////////////////////////////////////////////////////////////////////////// // int Gragg_Bulirsch_Stoer( double (*f)(double, double), double y[], // // double x, double h, double *h_next, double epsilon ) // // // // Description: // // This function solves the differential equation y'=f(x,y) with the // // initial condition y=y[0] at x. The value returned is y[1] which // // is the value of y evaluated at x + h, h_next is the step size used so // // accuracy is maintained. // // // // Arguments: // // double *f Pointer to the function which returns the slope at (x,y) of // // integral curve of the differential equation y' = f(x,y) // // which passes through the point (x0,y0) corresponding to the // // initial condition y(x0) = y0. // // double y[] On input y[0] is the initial value of y at x, and on output // // y[1] is the value of y at x + h. // // double x Initial value of x. // // double h Step size // // double *h_next Step size required to maintain the accuracy. // // double epsilon The tolerance. // // // // Return Values: // // The solution of y' = f(x,y) at x + h starting with y[0] at x. // // -1 failed to converge, -2 attempt to divide by zero . // //////////////////////////////////////////////////////////////////////////////// // // int Gragg_Bulirsch_Stoer( double (*f)(double, double), double y[], double x, double h, double *h_next, double epsilon ) { double step_size2[ATTEMPTS]; double aux_array[ATTEMPTS]; double dum; double est; double old_zero; int i; for (i = 0; i < ATTEMPTS; i++) { est = Graggs_Method( f, y[0], x, x+h, number_of_steps[i] ); step_size2[i] = (dum = h / (double) number_of_steps[i], dum * dum); if ( Rational_Function_Zero_Extrapolate( &y[1], step_size2, est, aux_array, i) < 0 ) return -2; if ( fabs(y[1] - old_zero) < epsilon ) { if ( i == LENGTH_SCALE ) *h_next = 0.95 * h; else if ( i == LENGTH_SCALE - 1 ) *h_next = 1.1 * h; else *h_next = (double) number_of_steps[LENGTH_SCALE-1] / (double) number_of_steps[ATTEMPTS - 1] * h; return 0; } old_zero = y[1]; } return -1; } //////////////////////////////////////////////////////////////////////////////// // // // Description: // // Gragg's method, also called the modified midpoint method, for solving // // a 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. // // //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// // 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 // // integral curve of the differential equation y' = f(x,y) // // which passes through the point (x0,y0) corresponding to the // // initial condition y(x0) = y0. // // double y0 The initial value of y at x0. // // double x0 Initial value of x. // // double x Final value of x. // // int number_of_steps The number_of_steps must be a positive even // // integer. The step size h is (x - x0) / number_of_steps. // // // // Return Values: // // The solution y(x) of the initial value problem y' = f(x,y) with // // y(x0) = y0. // // // //////////////////////////////////////////////////////////////////////////////// // // static 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) ); } //////////////////////////////////////////////////////////////////////////////// // int Rational_Function_Zero_Extrapolate( double *zero, double x[], // // double f, double aux_array[], int n ); // // // // Description: // // This function on f(x) and support points (x[i],f[x[i]]), this routine // // estimates the value of f(0) by collocating the // // // // The algorithm for extrapolating to x = 0 a sequence of support points // // (x[i],f[i]) for i = 0,...,n is to form a tableau T[row,col], // // row = 0,...,n and col = -1,0,...,n in which // // T[row,col] is defined inductively as follows: // // T[row,-1] = 0, T[row,0] = f[row], and // // T[row,col] = T[row,col-1] + // // ( (T[row,col-1] - T[row-1,col-1] )* (T[row,col-1] - T[row-1,col-2]) ) // // / ( ( x[row-col] / x[row] ) * (T[row-1,col-1] - T[row-1,col-2] ) // // - ( T[row,col-1] - T[row-1,col-2] ) ) // // // // Arguments: // // double *zero The estimation of f(0). // // double x[] The support abscissa. // // double f The support ordinate f[n]. // // double aux_array An auxiliary array containing the last row of the // // tableau described above. The dimension of aux_array // // must be the maximum number of support points. // // int n The index of the support ordinate. // // // // Return Values: // // If successful, the value returned is 0. If an attempt to divide by // // zero, -1 is returned. // // // //////////////////////////////////////////////////////////////////////////////// // // static int Rational_Function_Zero_Extrapolate( double *zero, double x[], double f, double aux_array[], int n ) { double back_two_columns; // T[row,col-2]; double old_aux; // T[row-1,col]; double new_value; // T[row,col]; double vertical_diff; // T[row,col]-T[row-1,col] double backslant_diff; // T[row,col]-T[row,col-1] double forwardslant_diff; // T[row,col]-T[row-1,col-1]; double denominator; int i; if (n == 0) { aux_array[0] = f; return 0; } if ( x[n] == 0.0 ) { *zero = f; return 0; } back_two_columns = 0.0; old_aux = aux_array[0]; aux_array[0] = f; for (i = 0; i < n; i++) { vertical_diff = aux_array[i] - old_aux; backslant_diff = aux_array[i] - back_two_columns; forwardslant_diff = backslant_diff - vertical_diff; denominator = x[n-i-1]/x[n] * forwardslant_diff - backslant_diff; if (denominator == 0.0) return -1; back_two_columns = old_aux; old_aux = aux_array[i+1]; aux_array[i+1] = aux_array[i] + vertical_diff * backslant_diff / denominator; } *zero = aux_array[n]; return 0; }