Newton's 3/8 Rule

Description

Newton's 3/8 rule approximates the integral of a function f(x) on the closed and bounded interval [a,a+h] of length h > 0 by the integral on [a,a+h] of the cubic passing through the points (a,f(a)), (a+h/3,f(a+h/3)), (a+2h/3,f(a+2h/3)) and (a+h,f(a+h)). The composite Newton's 3/8 rule is used to approximate the integral of a function f(x) over a closed and bounded interval [a,b] where a < b, by decomposing the interval [a,b] into n > 1 subintervals of equal length h = (b - a) / n, then adding the results of applying the Newton's 3/8 rule to each subinterval. By abuse of language both the composite Newton's 3/8 rule and Newton's 3/8 rule are referred to simply as Newton's 3/8 rule.
Let I_[a,b](f) be the integral of f(x) over the closed and bounded interval [a,b], and let N_h(f) be the result of applying the Newton's 3/8 rule with n subintervals of length h, i.e.
N_h(f)=(h/8) [ f(a) + 3f(a+h/3) + 3f(a+2h/3) + 2f(a+h) + ··· + 2f(b-h) + 3f(b-2h/3) + 3f(b-h/3) + f(b) ]. An immediate consequence of the Euler-Maclaurin summation formula relates I_[a,b](f) and N_h(f)
N_h(f) = I_[a,b](f) + (h⁴/6480) [ f'''(b) - f'''(a) ] - (h⁶/244944) [ f⁽⁵⁾(b) - f⁽⁵⁾(a) ] + ··· + K h^2p-2 [f^(2p-3)(b) - f^(2p-3)(a) ] + O(h^2p) ,where f''', f⁽⁵⁾, and f^(2p-3) are the third, fifth and (p-3)rd derivatives of f and K is a constant.

The last term, O(h^2p) is important. Given an infinitely differentiable function in which the first 2p-3 derivatives vanish at both endpoints of the interval of integration, it is not true that N_h(f) = I_[a,b](f) but rather what the theorem says is that
lim_h->0( N_h(f) - I_[a,b](f) ) / h^2p < M, where M > 0.

If f is at least four times differentiable on the interval [a,b], then applying the mean-value theorem to
N_h(f) - I_[a,b](f) = (h⁴/6480) [ f'''(b) - f'''(a) ] - (h⁶/244944) [ f⁽⁵⁾(b) - f⁽⁵⁾(a) ] + ··· + K h^2p-2 [f^(2p-3)(b) - f^(2p-3)(a) ] + O(h^2p) yields the standard truncation error expression
N_h(f) - I_[a,b](f) = (h⁴/6480) (b-a) f⁽⁴⁾(c), for some point c where a <= c <= b.A corollary of which is that if f⁽⁴⁾(x) = 0 for all x in [a,b], i.e. if f(x) is a cubic, then Newton's 3/8 rule is exact.

The Euler-Maclaurin summation formula also shows that usually n should be chosen large enough so that h = (b - a) / n < 1. For example, if h = 0.1 then
N_0.1(f) = I_[a,b](f) + 1.5·10^-8 [ f'''(b) - f'''(a) ] - 4.1·10^-12 [ f⁽⁵⁾(b) - f⁽⁵⁾(a) ] + ··· and if h = 0.01 then
N_0.01(f) = I_[a,b](f) + 1.5·10^-12 [ f'''(b) - f'''(a) ] - 4.1·10^-18 [ f⁽⁵⁾(b) - f⁽⁵⁾(a) ] + ··· while if h = 10 then
N₁₀(f) = I_[a,b](f) + 1.54 [ f'''(b) - f'''(a) ] - 4.08 [ f⁽⁵⁾(b) - f⁽⁵⁾(a) ] +···However, if the function f(x) is a cubic, then n may be chosen to be 1.

Besides the truncation error described above, the algorithm is also subject to the usual round-off errors and errors due to number of significant bits in the machine representation of floating point numbers.

The source code below is programmed in double precision, but the number of significant bits can easily be extended to extended precision by changing double to long double and by affixing an L to any floating point number e.g. 0.5 to 0.5L.

Mathematically addition is associative, (a+b)+c = a + (b+c), but if intermediate results are rounded to machine precision, machine addition is not necessarily associative. E.g. consider a decimal machine with 3 significant digits, then ( ( 0.400 + 0.400 ) + 0.400 ) + 122 = 123 > 122 = 0.400 + ( 0.400 + ( 0.400 + 122 ) ). In order to reduce the effect of intermediate result round-off errors when adding a series of floating point numbers, there are two versions of the rectangular rule, one which adds left to right N_h(f)=h/8 [ (((···( ( f(a) + 3f(a+h/3) ) + 3f(a+2h/3) ) + ··· + 3f(b-2h/3) ) + f(b)]and the other which adds right to leftN_h(f)=h/8 [ f(a) + ( 3f(a+h/3) + ( 3f(a+2h/3) + ··· + (f(b-2h/3) + f(b) )···))) ].In general, if the function is increasing in the interval of integration, addition should be performed left to right and if the function is decreasing in the interval of integration, addition should be performed right to left.

Function List

Functions:

double Newtons_3_8_Rule_Sum_LR( double a, double h, int n, double (*f)(double) )

Integrate the user supplied function (*f)(x) from a to a + nh where a is the lower limit of integration, h > 0 is the length of each subinterval, and n > 0 is the number of subintervals. The sum is performed from left to right.
double Newtons_3_8_Rule_Sum_RL( double a, double h, int n, double (*f)(double) )

Integrate the user supplied function (*f)(x) from a to a + nh where a is the lower limit of integration, h > 0 is the length of each subinterval, and n > 0 is the number of subintervals. The sum is performed from right to left.
double Newtons_3_8_Rule_Tab_Sum_LR( double h, int n, double f[] )

Integrate the function f[] given as an array of dimension n whose i^th element is the function evaluated at the midpoint of the i^th subinterval, where h > 0 is the length of each subinterval, and n > 0 is the number of subintervals. The sum is performed from left to right.
double Newtons_3_8_Rule_Tab_Sum_RL( double h, int n, double f[] )

Integrate the function f[] given as an array of dimension n whose i^th element is the function evaluated at the midpoint of the i^th subinterval, where h > 0 is the length of each subinterval h, and n > 0 is the number of subintervals. The sum is performed from right to left.

C Source

File: newtons_3_8_rule.c

////////////////////////////////////////////////////////////////////////////////
// File: newtons_3_8_rule.c                                                   //
// Routines:                                                                  //
//    Newtons_3_8_Rule_Sum_LR                                                 //
//    Newtons_3_8_Rule_Sum_RL                                                 //
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
//  Description:                                                              //
//     A (closed) Newton-Cotes formula is an interpolatory quadrature formula //
//     using equally spaced nodes for approximating the integral of a         //
//     function over a closed and bounded interval.                           //
//                                                                            //
//     Newton's 3/8ths rule is a (closed) Newton-Cotes formula for which the  //
//     interpolating polynomial is a cubic, the cubic passing through         //
//     (x,f(x)) at the endpoints of the interval and through (x1,f(x1)) and   //
//     (x2,f(x2)) where x1 is 1/3 of the distance from to left endpoint to    //
//     the right endpoint and x2 is 1/3 of the distance from the right enpoint//
//     to the left endpoint.                                                  //
//     The estimation of the integral then is simply the area beneath the     //
//     cubic from the left endpoint to the right endpoint.                    //
//                                                                            //
//     In order to integrate a function using Newton's 3/8ths rule 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 the Newton's 3/8ths rule to each      //
//     subinterval.  By abuse of language, this composite Newton's 3/8ths     //
//     rule is simply also called Newton's 3/8ths rule.                       //
//                                                                            //
//     Let I(f) be the approximation of the integral of f(x) over [a,b] then  //
//     Newton's 3/8ths rule is:                                               //
//       I(f) = h/8 { f(a) + 3f(a+h/3) + 3f(a+2h/3) + 2f(a+h) + ... +         //
//                                          3f(b-2h/3) + 3 f(b-h/3) + f(b) }  //
//     where h is the subinterval length, h = (b - a) / n, n is the number of //
//     subintervals.                                                          //
//                                                                            //
//     The truncation error estimate for Newton's 3/8ths rule is              //
//                    -  * n * h^5  f''''(xm) / 6480                          //
//     where f''''(xm) is the value of the fourth derivative evaluated at     //
//     some (unknown) point a <= xm <= b.                                     //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
//  double Newtons_3_8_Rule_Sum_LR( double a, double h, int n,                //
//                                                    double (*f)(double) );  //
//                                                                            //
//  Description:                                                              //
//     This function integrates f(x) from a to a+nh using the Newton's 3/8s   //
//     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, h > 0.                      //
//     int    n   The number of subintervals, the upper limit of integration  //
//                is a + n * h.                                               //
//     double *f  Pointer to the integrand, a function of a single variable   //
//                of type double.                                             //
//                                                                            //
//  Return Values:                                                            //
//     The integral of f(x) from a to (a + n * h).                            //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
double Newtons_3_8_Rule_Sum_LR( double a, double h, int n,
                                                       double (*f)(double) ) { 
   double h3 =  0.333333333333333333333 * h;
   double h23 = h3 + h3;
   double bound = a + ( n - 0.1666666666666666666666 ) * h;
   double integral = (*f)(a) + 3.0 * ( (*f)(a + h3) + (*f)(a + h23) );

   for (a += h; a < bound; a += h) 
      integral += 2.0 * (*f)(a) + 3.0 * ( (*f)(a + h3) + (*f)(a + h23) );
 
   return 0.125 * h * ( integral + (*f)(a) );
}


////////////////////////////////////////////////////////////////////////////////
//  double Newtons_3_8_Rule_Sum_RL( double a, double h, int n,          //
//                                                    double (*f)(double) );  //
//                                                                            //
//  Description:                                                              //
//     This function integrates f(x) from b-nh to b using the Newton's 3/8s   //
//     rule by summing from the right end of the interval to the left end.    //
//                                                                            //
//  Arguments:                                                                //
//     double a   The lower limit of the integration interval.                //
//     double h   The length of each subinterval, h > 0.                      //
//     int    n   The number of subintervals, the upper limit of integration  //
//                is a + n * h.                                               //
//     double *f  Pointer to the integrand, a function of a single variable   //
//                of type double.                                             //
//                                                                            //
//  Return Values:                                                            //
//     The integral of f(x) from a to (a + n * h).                            //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
double Newtons_3_8_Rule_Sum_RL( double a, double h, int n,
                                                       double (*f)(double) ) { 
   double h3 =  0.333333333333333333333 * h;
   double h23 = h3 + h3;
   double bound = a + 0.5 * h3;
   double x = a + n * h;
   double integral = (*f)(x) + 3.0 * ( (*f)(x - h3) + (*f)(x - h23) );

   for (x -= h; x > bound; x -= h) 
      integral += 2.0 * (*f)(x) + 3.0 * ( (*f)(x - h3) + (*f)(x - h23) );
 
   return 0.125 * h * ( integral + (*f)(a) );
}

File: simpsons_rule_tab.c

////////////////////////////////////////////////////////////////////////////////
// File: newtons_3_8_rule_tab.c                                               //
// Routines:                                                                  //
//    Newtons_3_8_Rule_Tab_Sum_LR                                             //
//    Newtons_3_8_Rule_Tab_Sum_RL                                             //
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
//  Description:                                                              //
//     A (closed) Newton-Cotes formula is an interpolatory quadrature formula //
//     using equally spaced nodes for approximating the integral of a         //
//     function over a closed and bounded interval.                           //
//                                                                            //
//     Newton's 3/8ths rule is a (closed) Newton-Cotes formula for which the  //
//     interpolating polynomial is a cubic, the cubic passing through         //
//     (x,f(x)) at the endpoints of the interval and through (x1,f(x1)) and   //
//     (x2,f(x2)) where x1 is 1/3 of the distance from to left endpoint to    //
//     the right endpoint and x2 is 1/3 of the distance from the right enpoint//
//     to the left endpoint.                                                  //
//     The estimation of the integral then is simply the area beneath the     //
//     cubic from the left endpoint to the right endpoint.                    //
//                                                                            //
//     In order to integrate a function using Newton's 3/8ths rule 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 the Newton's 3/8ths rule to each      //
//     subinterval.  By abuse of language, this composite Newton's 3/8ths     //
//     rule is simply also called Newton's 3/8ths rule.                       //
//                                                                            //
//     Let I(f) be the approximation of the integral of f(x) over [a,b] then  //
//     Newton's 3/8ths rule is:                                               //
//       I(f) = h/8 { f(a) + 3f(a+h/3) + 3f(a+2h/3) + 2f(a+h) + ... +         //
//                                          3f(b-2h/3) + 3 f(b-h/3) + f(b) }  //
//     where h is the subinterval length, h = (b - a) / n, n is the number of //
//     subintervals.                                                          //
//                                                                            //
//     The truncation error estimate for Newton's 3/8ths rule is              //
//                    -  * n * h^5  f''''(xm) / 6480                          //
//     where f''''(xm) is the value of the fourth derivative evaluated at     //
//     some (unknown) point a <= xm <= b.                                     //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
//  double Newtons_3_8_RuleTab_Sum_LR( double h, int n, double f[] );         //
//                                                                            //
//  Description:                                                              //
//     Assuming that a is the lower limit of integration and the ith element  //
//     of the array f[] is the value of the function evaluated at the left    //
//     endpoint of the ith subinterval for 0 <= i < n and f[n] is the right   //
//     endpoint of the nth interval, i.e. f[i] = f( a + i*h ), i = 0,...,n,   //
//     then  this routine integrates f(x) using Newton's 3/8 rule by summing  //
//     from the left end of the interval to the right end.                    //
//                                                                            //
//  Arguments:                                                                //
//     double h   The length of each subinterval, h > 0.                      //
//     int    n   The number of subintervals, the upper limit of integration  //
//                is a + n * h, where a is the lower limit of integration.    //
//     double f[] Array of function values the ith element of which is the    //
//                f[i] = f(a + i*h/3), i = 0,...,3n, i.e. the dimension of f  //
//                is 3n + 1.                                                  //
//                                                                            //
//  Return Values:                                                            //
//     The integral of f(x) from a to (a + n * h).                            //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
double Newtons_3_8_Rule_Tab_Sum_LR( double h, int n, double f[] ) {
 
   double integral = f[0] + 3.0 * ( f[1] + f[2] );
   int i;
   
   for (i = 3; i < ( n + n + n ); i += 3)
      integral += 2.0 * f[i] + 3.0 * ( f[i+1] + f[i+2] );

   return 0.125 * h * ( integral + f[i] );
}


////////////////////////////////////////////////////////////////////////////////
//  double Newtons_3_8_RuleTab_Sum_RL( double h, int n, double f[] );         //
//                                                                            //
//  Description:                                                              //
//     Assuming that a is the lower limit of integration and the ith element  //
//     of the array f[] is the value of the function evaluated at the left    //
//     endpoint of the ith subinterval for 0 <= i < n and f[n] is the right   //
//     endpoint of the nth interval, i.e. f[i] = f( a + i*h ), i = 0,...,n,   //
//     then  this routine integrates f(x) using Newton's 3/8 rule by summing  //
//     from the right end of the interval to the left end.                    //
//                                                                            //
//  Arguments:                                                                //
//     double h   The length of each subinterval, h > 0.                      //
//     int    n   The number of subintervals, the upper limit of integration  //
//                is a + n * h, where a is the lower limit of integration.    //
//     double f[] Array of function values the ith element of which is the    //
//                f[i] = f(a + i*h/3), i = 0,...,3n, i.e. the dimension of f  //
//                is 3n + 1.                                                  //
//                                                                            //
//  Return Values:                                                            //
//     The integral of f(x) from a to (a + n * h).                            //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
double Newtons_3_8_Rule_Tab_Sum_RL( double h, int n, double f[] ) { 
  
   int i = n + n + n; 
   double integral = f[i] + 3.0 * ( f[i-1] + f[i-2] );

   for (i -= 3 ; i > 0; i -= 3) 
      integral += 2.0 * f[i] + 3.0 * ( f[i-1] + f[i-2] );

   return 0.125 * h * ( integral + f[0] );
}