Hermite Quadrature
C Source
////////////////////////////////////////////////////////////////////////////////
// File: hermite_quadrature_1_derivative.c //
// Routines: //
// Hermite_Quadrature_1_Derivative_LR //
// Hermite_Quadrature_1_Derivative_RL //
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// //
// Description: //
// The Hermite quadrature formula extends the trapezoidal rule by using //
// the values of the first derivatives of the integrand at the endpoints //
// of the integration interval. This technique is convenient when the //
// derivatives of the integrand are easy to calculate. //
// //
// This version of the Hermite quadrature is sometimes called the //
// corrected trapezoidal rule. //
// //
// The integral over a fundamental interval [0,h] is //
// (h/2) [f(0) + f(h) + (h/6) (f'(0) - f'(h) ] //
// with a local truncation error of //
// h^5 / 720 f''''(xm). //
// //
// In order to integrate a function using Hermite quadrature over a //
// closed and bounded interval [a,b], divide the interval [a,b] into N //
// subintervals of length h = (b-a) / N. The integral of the function //
// over the interval [a,b] is the sum of the integrals of the function //
// over the subintervals [a,a+h], [a+h, a+2h],...,[b-h,b]. The integral //
// is approximated then by applying Hermite quadrature to each //
// subinterval. //
// //
// Let I(f) be the approximation of the integral of f(x) over [a,b] then //
// Hermite quadrature is: //
// I(f) = h { f(x[0]) / 2 + f(x[1]) + ... + f(x[n-1]) + f(x[n]) / 2) } //
// + h^2/12 (f'(a) - f'(b)) //
// where h is the step size, h = x[i] - x[i-1]; n is the number of //
// steps, n = abs(b - a) / h and x[0] = a, x[i] = x[i-1] + h for //
// i = 1,...,n-1 so that x[n] = b. //
// //
// The truncation error estimate for the Hermite quadrature is //
// n * h^5 f''''(xm) / 720 //
// where f''''(xm) is the value of the fourth derivative evaluated at //
// some (unknown) point a <= xm <= b. //
// //
// Note! If f'(b) = f'(a), Hermite quadrature is equivalent to the //
// trapezoidal rule. In this case, the trapezoidal rule is superior to //
// Simpson's rule. //
// //
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// double Hermite_Quadrature_1_Derivative_LR( double a, double h, int n, //
// double (*f)(double), double (*df)(double) ); //
// //
// Description: //
// This function integrates f(x) from a to a+nh using the corrected //
// trapezoidal rule by summing from the left end of the interval to the //
// right end. //
// //
// Arguments: //
// double a The lower limit of the integration interval. //
// double h The length of each subinterval. //
// int n The number of steps, the upper limit of integration is //
// a + n * h. //
// double *f Pointer to the integrand, a function of a single variable //
// of type double. //
// double *df Pointer to the derivative of the integrand, a function of a //
// single variable of type double. //
// //
// Return Values: //
// The integral of f(x) from a to (a + n * h). //
// //
// Example: //
////////////////////////////////////////////////////////////////////////////////
// //
double Hermite_Quadrature_1_Derivative_LR( double a, double h, int n,
double (*f)(double), double (*df)(double) ) {
double upper_limit = a + (n - .5) * h;
double integral = 0.5 * (*f)(a)
+ 0.0833333333333333333333 * h * ( (*df)(a) - (*df)(a+n*h) );
for (a += h; a < upper_limit; a += h) integral += (*f)(a);
integral += 0.5 * (*f)(a);
return h * integral;
}
////////////////////////////////////////////////////////////////////////////////
// double Hermite_Quadrature_1_Derivative_RL( double b, double h, int n, //
// double (*f)(double), double (*df)(double) ); //
// //
// Description: //
// This function integrates f(x) from b-nh to b using the corrected //
// trapezoidal rule by summing from the left end of the interval to the //
// right end. //
// //
// Arguments: //
// double b The upper limit of the integration interval. //
// double h The length of each subinterval. //
// int n The number of steps, the lower limit of integration is //
// b - n * h. //
// double *f Pointer to the integrand, a function of a single variable //
// of type double. //
// double *df Pointer to the derivative of the integrand, a function of a //
// single variable of type double. //
// //
// Return Values: //
// The integral of f(x) from (b - n * h) to b. //
// //
// Example: //
////////////////////////////////////////////////////////////////////////////////
// //
double Hermite_Quadrature_1_Derivative_RL( double b, double h, int n,
double (*f)(double), double (*df)(double) ) {
double lower_limit = b - (n - .5) * h;
double integral = 0.5 * (*f)(b)
+ 0.0833333333333333333333 * h * ( (*df)(b-n*h) - (*df)(b) );
for (b -= h; b > lower_limit; b -= h) integral += (*f)(b);
integral += 0.5 * (*f)(b);
return h * integral;
}