Ralston's Fourth Order Method
Description
Ralston's fourth order method is a Runge-Kutta method for approximating the solution of the initial value problem y'(x) = f(x,y); y(x0) = y0 which evaluates the integrand,f(x,y), four times per step. See the comments in the source code for the algorithm.
This method is a fourth order procedure for which Richardson extrapolation can be used.
Function List
Functions:
- double Runge_Kutta_Ralston_4_Method( double (*f)(double, double), double y0, double x0, double h, int number_of_steps )
This function uses Ralston's fourth order method to return the estimate of the solution of the initial value problem, y' = f(x,y); y(x0) = y0, at x = x0 + h * n, where h is the step size and n is number_of_steps.
- double Runge_Kutta_Ralston_4_Richardson( double (*f)(double, double), double y0,
double x0, double h, int number_of_steps, int richardson_columns )
This function uses Ralston's fourth order method together with Richardson extrapolation to the limit to return the estimate of the solution of the initial value problem, y' = f(x,y); y(x0) = y0, at x = x0 + h * n, where h is the step size and n is number_of_steps. The argument richardson_columns is
the number of step size halving + 1 used in Richardson extrapolation so that if richardson_columns = 1 then no extrapolation to the limit is performed.
- void Runge_Kutta_Ralston_4_Integral_Curve( double (*f)(double, double), double y[ ], double x0, double h, int number_of_steps_per_interval, int number_of_intervals )
This function uses Ralston's fourth order method to estimate the solution of the initial value problem, y' = f(x,y); y(x0) = y0, at x = x0 + h * n * m, where h is the step size and n is the interval number 0 <= n <= number_of_intervals, and m is the number_of_steps_per_interval. The values are return in the array y[ ] i.e. y[n] = y(x0 + h * m * n), where m, n are as above.
- void Runge_Kutta_Ralston_4_Richardson_Integral_Curve( double (*f)(double, double), double y[ ], double x0, double h, int number_of_steps_per_interval, int number_of_intervals, int richardson_columns )
This function uses Ralston's fourth order method together with Richardson extrapolation to the limit to estimate the solution of the initial value problem, y' = f(x,y); y(x0) = y0, at x = x0 + h * n * m, where h is the step size and n is the interval number 0 <= n <= number_of_intervals, and m is the number_of_steps_per_interval. The values are return in the array y[ ] i.e. y[n] = y(x0 + h * m * n), where m, n are as above. The argument richardson_columns is
the number of step size halving + 1 used in Richardson extrapolation so that if richardson_columns = 1 then no extrapolation to the limit is performed.
C Source
////////////////////////////////////////////////////////////////////////////////
// File: runge_kutta_ralston_4.c //
// Routines: //
// Runge_Kutta_Ralston_4_Method //
// Runge_Kutta_Ralston_4_Richardson //
// Runge_Kutta_Ralston_4_Integral_Curve //
// Runge_Kutta_Ralston_4_Richardson_Integral_Curve //
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// //
// Description: //
// The fourth order Ralston's method is a four stage Runge-Kutta method //
// for approximating the solution of the differential equation //
// y'(x) = f(x,y) with initial condition y(x0) = c. The fourth order //
// Ralston's method evaluates f(x,y) four times per step. For step i+1, //
// y[i+1] = y[i] + (263 + 24s5)/1812 k1 + (125-1000s5)/3828 k2 //
// + 1024*(3346 + 1623s5)/5924787 k3 + (30 - 4s5)/123 k4, //
// where //
// k1 = h * f(x[i],y[i]), //
// k2 = h * f(x[i]+2/5 h, y[i]+2/5 k1), //
// k3 = h * f(x[i] + (7/8 - 3s5/16)h, y[i] + (-2889+1428s5)/1024 k1 //
// + (3785-1620s5)/1024 k2) //
// k4 = h * f(x[i] +h, y[i] + (-3365+2094s5)/6040 k1 //
// (-975-3046s5)/2552 k2 //
// (467040+203968s5)/2400845 k3) ) //
// and s5 = sqrt(5), x[i] = x0 + i * h. //
// //
// This version of the Runge-Kutta method is a fourth order method. //
// Richardson extrapolation can be used to increase the order and //
// accuracy. //
// //
////////////////////////////////////////////////////////////////////////////////
#define SQRT5 2.236067977499789694091736687312762
static const double a2 = 0.40;
static const double a3 = (14.0 - 3.0 * SQRT5) / 16.0;
static const double b31 = (-2889.0 + 1428.0 * SQRT5) / 1024.0;
static const double b32 = (3785.0 - 1620.0 * SQRT5) / 1024.0;
static const double b41 = (-3365.0 + 2094.0 * SQRT5) / 6040.0;
static const double b42 = (-975.0 - 3046.0 * SQRT5) / 2552.0;
static const double b43 = (467040.0 + 203968.0 * SQRT5) / 240845.0;
static const double g1 = (263.0 + 24.0 * SQRT5) / 1812.0;
static const double g2 = (125.0 - 1000.0 * SQRT5) / 3828.0;
static const double g3 = 1024.0 * (3346.0 + 1623.0 * SQRT5) / 5924787.0;
static const double g4 = (30.0 - 4.0 * SQRT5) / 123.0;
////////////////////////////////////////////////////////////////////////////////
// double Runge_Kutta_Ralston_4_Method( double (*f)(double, double), //
// double y0, double x0, double h, int number_of_steps ); //
// //
// Description: //
// This routine uses the 4th order Runge-Kutta method described above to //
// approximate the solution at x = x0 + h * number_of_steps of the initial//
// value problem y'=f(x,y), y(x0) = y0. //
// //
// 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 //
// passes through the point (x0,y0). //
// double y0 //
// The initial value of y at x = x0. //
// double x0 //
// The initial value of x. //
// double h //
// The step size. //
// int number_of_steps //
// The number of steps. Must be a nonnegative integer. //
// //
// Return Values: //
// The solution of the initial value problem y' = f(x,y), y(x0) = y0 at //
// x = x0 + number_of_steps * h. //
// //
////////////////////////////////////////////////////////////////////////////////
// //
double Runge_Kutta_Ralston_4_Method( double (*f)(double, double), double y0,
double x0, double h, int number_of_steps ) {
double k1, k2, k3, k4;
double h2 = a2 * h;
double h3 = a3 * h;
while ( --number_of_steps >= 0 ) {
k1 = (*f)(x0,y0);
k2 = (*f)(x0+h2, y0 + h2 * k1);
k3 = (*f)(x0+h3, y0 + h * (b31 * k1 + b32 * k2));
x0 += h;
k4 = (*f)(x0, y0 + h * (b41 * k1 + b42 * k2 + b43 * k3));
y0 += h * ( g1 * k1 + g2 * k2 + g3 * k3 + g4 * k4 );
}
return y0;
}
////////////////////////////////////////////////////////////////////////////////
// double Runge_Kutta_Ralston_4_Richardson( double (*f)(double, double), //
// double y0, double x0, double h, int number_of_steps, //
// int richardson_columns) //
// //
// Description: //
// This routine uses the 4th order Runge-Kutta method described above //
// together with Richardson extrapolation to approximate the solution at //
// x = x0 + h * number_of_steps of the initial value problem y'=f(x,y), //
// y(x0) = y0. //
// //
// 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 //
// passes through the point (x0,y0). //
// double y0 //
// The initial value of y at x = x0. //
// double x0 //
// The initial value of x. //
// double h //
// The step size. //
// int number_of_steps //
// The number of steps. Must be nonnegative. //
// int richardson_columns //
// The maximum number of columns to use in the Richardson //
// extrapolation to the limit. //
// //
// Return Values: //
// The solution of the initial value problem y' = f(x,y), y(x0) = y0 at //
// x = x0 + number_of_steps * h. //
// //
////////////////////////////////////////////////////////////////////////////////
// //
static const double richardson[] = {
1.0 / 15.0, 1.0 / 31.0, 1.0 / 63.0, 1.0 / 127.0, 1.0 / 255.0, 1.0 / 511.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 Runge_Kutta_Ralston_4_Richardson( double (*f)(double, double), double y0,
double x0, double h, int number_of_steps, int richardson_columns ) {
double dt[MAX_COLUMNS]; // dt[i] is the last element in column i.
double integral, delta, h_used;
int j,k, number_sub_intervals;
richardson_columns = max(1, min(MAX_COLUMNS, richardson_columns));
while ( --number_of_steps >= 0 ) {
h_used = h;
number_sub_intervals = 1;
for (j = 0; j < richardson_columns; j++) {
integral = Runge_Kutta_Ralston_4_Method( f, y0, x0, h_used,
number_sub_intervals);
for ( k = 0; k < j; k++) {
delta = integral - dt[k];
dt[k] = integral;
integral += richardson[k] * delta;
}
dt[j] = integral;
h_used *= 0.5;
number_sub_intervals += number_sub_intervals;
}
y0 = integral;
x0 += h;
}
return y0;
}
////////////////////////////////////////////////////////////////////////////////
// void Runge_Kutta_Ralston_4_Integral_Curve( double (*f)(double, double), //
// double y[], double x0, double h, int number_of_steps_per_interval, //
// int number_of_intervals ); //
// //
// Description: //
// This routine uses the 4th order Runge-Kutta method described above to //
// approximate the solution of the differential equation y'=f(x,y) with //
// the initial condition y = y[0] at x = x0. The values are returned in //
// y[n] which is the value of y evaluated at x = x0 + n * m * h, where m //
// is the number of steps per interval and n is the interval number, //
// 0 <= n <= number_of_intervals. //
// //
// 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,y[0]). //
// double y[] //
// On input y[0] is the initial value of y at x = x0. On output //
// for n >= 1, y[n] is the appproximation of the solution y(x) of //
// the initial value problem y' = f(x,y), y(x0) = y[0], at //
// x = x0 + nmh, where m is the number of steps per interval and n //
// is the interval number. //
// double x0 //
// Initial value of x. //
// double h //
// Step size //
// int number_of_steps_per_interval //
// The number of steps of length h used to calculate y[i+1] //
// starting from y[i]. //
// int number_of_intervals //
// The number of intervals, y[] should be dimensioned at least //
// number_of_intervals + 1. //
// //
// Return Values: //
// This routine is of type void and hence does not return a value. //
// The solution of y' = f(x,y) from x = x0 to x = x0 + n * m * h, //
// where n is the number of intervals and m is the number of steps per //
// interval, is stored in the input array y[]. //
// //
////////////////////////////////////////////////////////////////////////////////
// //
void Runge_Kutta_Ralston_4_Integral_Curve( double (*f)(double, double),
double y[], double x0, double h, int number_of_steps_per_interval,
int number_of_intervals ) {
double k1, k2, k3, k4;
double h2 = a2 * h;
double h3 = a3 * h;
int i;
while ( --number_of_intervals >= 0 ) {
*(y+1) = *y;
y++;
for (i = 0; i < number_of_steps_per_interval; i++) {
k1 = (*f)(x0,*y);
k2 = (*f)(x0+h2, *y + h2 * k1);
k3 = (*f)(x0+h3, *y + h * (b31 * k1 + b32 * k2));
x0 += h;
k4 = (*f)(x0, *y + h * (b41 * k1 + b42 * k2 + b43 * k3));
*y += h * ( g1 * k1 + g2 * k2 + g3 * k3 + g4 * k4 );
}
}
}
////////////////////////////////////////////////////////////////////////////////
// void Runge_Kutta_Ralston_4_Richardson_Integral_Curve( //
// double (*f)(double, double), double y[], double x0, double h, //
// int number_of_steps_per_interval, int number_of_intervals, //
// int richardson_columns ) //
// //
// Description: //
// This routine uses the 4th order Runge-Kutta method described above //
// together with Richardson extrapolation to the limit (as h -> 0) to //
// approximate the solution of the differential equation y'=f(x,y) with //
// the initial condition y = y[0] at x = x0. //
// The values are returned in y[], in which y[n] is the value of y //
// evaluated at x = x0 + n * m * h, where m is the number of steps per //
// interval and n is the interval number, 0 <= n <= number_of_intervals. //
// //
// 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,y[0]). //
// double y[] //
// On input y[0] is the initial value of y at x = x0. On output //
// for n >= 1, y[n] is the appproximation of the solution y(x) of //
// the initial value problem y' = f(x,y), y(x0) = y[0], at //
// x = x0 + nmh, where m is the number of steps per interval and n //
// is the interval number. //
// double x0 //
// Initial value of x. //
// double h //
// Step size //
// int number_of_steps_per_interval //
// The number of steps of length h used to calculate y[i+1] //
// starting from y[i]. //
// int number_of_intervals //
// The number of intervals, y[] should be dimensioned at least //
// number_of_intervals + 1. //
// int richardson_columns //
// The maximum number of columns to use in the Richardson //
// extrapolation to the limit. //
// //
// Return Values: //
// This routine is of type void and hence does not return a value. //
// The solution of y' = f(x,y) from x = x0 to x = x0 + n * m * h, //
// where n is the number of intervals and m is the number of steps per //
// interval, is stored in the input array y[]. //
// //
////////////////////////////////////////////////////////////////////////////////
// //
void Runge_Kutta_Ralston_4_Richardson_Integral_Curve(
double (*f)(double, double), double y[], double x0, double h,
int number_of_steps_per_interval, int number_of_intervals,
int richardson_columns ) {
double mh = (double) number_of_steps_per_interval * h;
while ( --number_of_intervals >= 0 ) {
*(++y) = Runge_Kutta_Ralston_4_Richardson( f, *y, x0, h,
number_of_steps_per_interval, richardson_columns );
x0 += mh;
}
return;
}