Simpson’s 1/3rd Rule

A brief introduction to the Simpson’s 1/3^rd rule and a uniform interval Composite Simpson’s 1/3^rd Rule implementation.

Simpson’s 1/3^rd Rule

The Simpson’s 1/3 rule is a numerical method to find the integral within some finite limits and . Simpson’s 1/3^rd rule approximates with a polynomial of degree two , ie. a parabola between the two limits and , and then finds the integral of that bounded parabola, and is used to represent the approximate integral . The integral of the approximated function is the area under the parabola bounded by the points and and by the positive side of the x axis. The quadratic function has three points common to the function , which are as follows: The end points of the approximate quadratic function is the same as the function at , . And takes the same value of the function at point . Thus three points are fixed each in equal interval and a parabola is drawn through these three points and the area under the parabola through these points bounded by and and the positive side of the X axis is found, which is used as the approximated integral value. The below iterative formula an be used to find the integral of a function using Simpson’s 1/3^rd rule.

The above formula approximates with only one parabola. To represent the function more accurately with with method, the function is divided in many intervals say having same width say , such that and then a parabola is fit within each of the intervals to approximate the integrals of the function within those intervals. The evaluation of the integral within each such interval using the above single interval iterative formula is used to find the overall evaluation of the integral within the limits and . This is known as the Simpson’s 1/3^rd Composite Rule. The iterative method for the Simpson’s 1/3^rd Composite Rule can be found by repetitively applying the single interval formula in each intervals, and is found to be:

The Process

The Simpson’s Composite Rule function will calculate the first and the last ordinate separately and also calculate and sum up the odd and the even ordinates separately for better clarity, and at last the Simpson’s Composite Rule iterative formula would be used directly.
For demonstration purpose some sample functions are used:

1. 
2. 

FUNCTION simpson_13(f, a, b, n)
  Find width of each interval h = (b-a)/n 
  Find y0 = f(a), yn = f(b)
  Initialize x = a+h, i = 1, sum = 0
  Initialize even = 0, odd = 0
  /* Calculate the sum of the even and odd ordinates from (x + h) to (x + (n-1)h) */

  WHILE (i < n)
   DO
    IF i is even
     THEN
      even = even + f (x)
    ELSE IF i is odd
     THEN
      odd = odd + f(x)
    ENDIF
    x = x + h
  ENDWHILE

  /* Apply Simpson's 1/3rd Composite Formula to find the integral */

  Is = (h/3* {y0 + yn + 2 * even + 4 * odd})

ENDFUNCTION

Sourcecode

#include <stdio.h>
#include <math.h>

double f1 (double x);
double f2 (double x);

double simpson_13 (double (*f) (double x), double a, double b, int n);

int
main (void)
{
  double a, b, n;
  double Is;

  printf ("\nEnter a,b,n: ");
  scanf ("%lf %lf %lf", &a, &b, &n);

  printf ("\nf(x) = sin (2x) / (1+x)^5");

  /* Show integral computed with simpson's 1/3 rd rule */
  Is = simpson_13 (f1, a, b, n);
  printf ("\n\tI_Simpson_13rd (f(x), %g, %g, %g) = %g", a, b, n, Is);

  printf ("\n");

  printf ("\nf(x) = (1/x) + 5 + 10x^2;");

  /* Show integral computed with simpson's 1/3 rd rule */
  Is = simpson_13 (f2, a, b, n);
  printf ("\n\tI_Simpson_13rd (f(x), %g, %g, %g) = %g", a, b, n, Is);

  printf ("\n");

  return 0;
}

/* Takes a function pointer which should represent the single variable
 * function of which the integral is to be found. The function should
 * receive a double and return a double. The limits a and b and the
 * number of intervals n is provided. The value returned is the integral
 * of the passed function is the argument
 */
double
simpson_13 (double (*f) (double x), double a, double b, int n)
{
  double h;
  double y = 0, x, even = 0, odd = 0, y0, yn;
  int i;

  h = (b - a) / n;

  y0 = f (a);
  yn = f (b);

  for (i = 1, x = a + h; i < n; x = x + h, i++)
    {
      if (i % 2 == 0)
        even = even + f (x);
      else
        odd = odd + f (x);
    }

  y = (h / 3) * ((y0 + yn) + 2 * even + 4 * odd);

  return y;
}

double
f1 (double x)
{
  return sin (2 * x) / pow ((1 + x), 5);
}

double
f2 (double x)
{
  return (1 / x) + 5 + 10 * x * x;
}

Description

The two given function defined in the previous section is implemented with the functions f1 () and f2 (). These two functions accept a double type data, and returns double.

• double f1 (double x) : This function implements the equation sin (2 * x) / pow ((1 + x), 5) and returns the evaluated value.

• double f2 (double x) : This function implements the equation (1 / x) + 5 + 10 * x * x and returns the evaluated value.

• double simpson_13 (double (*f) (double x), double a, double b, int n) : This takes a pointer to a function pointer which accepts one double type data and returns a double type data. The next argument is double a, which is interpreted as the lower limit of the integral, the next argument double b is considered as the upper limit of the integral. int n is the number of intervals to be taken. After finding the width of each interval, this first calculates the first and the last ordinate y0 = f(a) and y1 = f(b). Then it enters the for loop and computes the values of the passed function starting from up to ). The values of the function values for the odd interval number are summed in a separate variable odd, and the values for the even interval numbers are summed in another variable even. At the end of the interval the composite iterative Simpson’s formula is directly used to find the integral, and the value is returned:

  (h / 3) * ((y0 + yn) + 2 * even + 4 * odd);

• int main (void) : The main function prompts the user to enter the upper and the lower limits and the number of intervals to be taken and calls the functions with the proper parameters. The simpson_13() is called as Is = simpson_13 (f1, a, b, n); which computes the integral of the function f1 within upper and lower limit a and b with the Simpson’s 1/3 rule.

Error

When is continuous and is the upper bound for the values
of on , then Simpson’s 1/3^rd is found to be:

The magnitude of the Simpson’s rule error decreases as the the step size is decreased as the function is better approximated with the second degree equation.

Comments

The Simpson’s 1/3^rd evaluation function is implemented to be independent of the function of which the integral is to be found, by allowing passing of the function to be integrated as a parameter. This implementation could be incorporated into an algebraic equation parser, which can be used to evaluate some user entered functions to be evaluated like for example:
ans = simpson_13 (parse_expression (user_expression, x), a, b, n) ; the ‘user_expression’ is being considered a single variable function with the variable replaced with the value of ‘x’ .

Output

Some sample outputs of the above code is shown below

Run 1: input = lower limit a = 1, upper limit b = 2, intervals n = 10

Enter a,b,n: 1 2 10

f(x) = sin (2x) / (1+x)^5
	I_Simpson_13rd (f(x), 1, 2, 10) = 0.00505831

f(x) = (1/x) + 5 + 10x^2;
	I_Simpson_13rd (f(x), 1, 2, 10) = 29.0265

Run 2: input = lower limit a = 5.3, upper limit b = 10.23, intervals n = 100

Enter a,b,n: 5.3 10.23 100

f(x) = sin (2x) / (1+x)^5
	I_Simpson_13rd (f(x), 5.3, 10.23, 100) = -3.20359e-05

f(x) = (1/x) + 5 + 10x^2;
	I_Simpson_13rd (f(x), 5.3, 10.23, 100) = 3097.71

Links

1. Check for more: http://en.wikipedia.org/wiki/Simpson%27s_rule
2. http://mathworld.wolfram.com/SimpsonsRule.html
3. Trapezoidal Rule: https://phoxis.org/2011/02/27/trapezoidal-rule/

References

1. Calculus and Analytic Geometry: Thomas Finney
2. Images from : Wikipedia