Locally the trapezoidal method is a third order method and therefore globally a second order method.
The trapezoidal method is a stable and convergent method with region of absolute stability
As a rule, the trapezoidal method is convenient if the function g(x,h,u) is simple.
////////////////////////////////////////////////////////////////////////////////
// File: trapezoidal_method.c //
// Routines: //
// Trapezoidal_Method //
// Trapezoidal_Integral_Curve //
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// //
// Description: //
// The trapezoidal method is an implicit method for approximating the //
// solution of the differential equation y'(x) = f(x,y) with initial //
// condition y(a) = c. The trapezoidal method can be derived by //
// expanding y(x+h) in a Taylor series about x, //
// y(x+h) = y(x) + h*y'(x) + (h^2 / 2)*y''(x) + (h^3 / 6)*y'''(x) +... //
// and substituting (y'(x+h)-y'(x))/h for y''(x) and y'(x) = f(x,y(x)), //
// y(x+h) = y(x) + (h/2) * (f(x+h,y(x+h) - f(x,h)) + O(h^3). //
// //
// Let y[n] = y(a + nh), x[n] = a + nh, then the recursion formula for //
// the trapezoidal method is: //
// y[n+1] = y[n] + (h/2) * ( f(x[n],y[n]) + f(x[n+1],y[n+1]) ) //
// //
// I.e. y[n+1] is the solution y to: //
// y - (h/2)*(f(x[n+1],y) = y[n] + (h/2)*f(x[n],y[n]). //
// //
// Thus locally the trapezoidal method is a third order method and //
// globally a second order method. //
// //
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// double Trapezoidal_Method( double (*f)(double, double), //
// double (*g)(double,double,double), double y0, double x0, //
// double h, int number_of_steps ); //
// //
// Description: //
// This routine uses the trapezoidal method 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 *g //
// Pointer to the function g(x,h,u) which returns the value y //
// such that y - h * f(x,y) = u. //
// 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 Trapezoidal_Method( double (*f)(double, double),
double (*g)(double,double,double), double y0, double x0,
double h, int number_of_steps ) {
double u;
double h2 = 0.5 * h;
while ( --number_of_steps >= 0 ) {
u = y0 + h2 * f(x0,y0);
x0 += h;
y0 = g(x0, h2, u);
}
return y0;
}
////////////////////////////////////////////////////////////////////////////////
// void Trapezoidal_Integral_Curve( double (*f)(double, double), //
// double (*g)(double,double,double), double y[], double x0, //
// double h, int number_of_steps_per_interval, //
// int number_of_intervals ); //
// //
// Description: //
// This routine uses the trapezoidal method 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 *g //
// Pointer to the function g(x,h,u) which returns the value y //
// such that y - h/2 f(x,y) = u. //
// double y[] //
// On input y[0] is the initial value of y at x = x0. On output //
// y[n+1] = y[n] + m*h*f(a + n*m*h, y[n] ), 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 Trapezoidal_Integral_Curve( double (*f)(double, double),
double (*g)(double,double,double), double y[], double x0, double h,
int number_of_steps_per_interval, int number_of_intervals ) {
double u;
double h2 = 0.5 * h;
int i;
while ( --number_of_intervals >= 0 ) {
*(y+1) = *y;
y++;
for (i = 0; i < number_of_steps_per_interval; i++) {
u = *y + h2 * f(x0,*y);
x0 += h;
*y = g(x0, h2, u);
}
}
}