COST (GBP)

0.00

Maths › Calculus › Quadrature ›

simpson

Computes the definite integral of a function using Simpson's first rule.

Controller: CodeCogs

Contents

Interface
Simpson
Simpson
Page Comments

Interface

C++

Overview

This module computes the area beneath either a user supplied function or a set of discrete points, using an approximation which assumes the function is shaped as a parabola between each two consecutive points.

{The background derivation of Simpon's first rule is given here:} Simpon's first rule represents the function between $\inline -h<=x<=h$ as a second order polynomial. The function is described by three ordinates, $\inline f_{-1}$ , $\inline f_0$ and $\inline f_1$ . In order to identify the truncation error, a fourth order polynomial is used,

$f(x) = f_0 + a x + b x^2 + c x^3 + d x ^4$

(2)

Once integrated this becomes,

$I = \int_{-h}^h f(x) \;dx = \left [ f_0 x + \frac{1}{2}a x^2 + \frac{1}{3} b x^3 + \frac{1}{4} c x^4 + \frac{1}{5} d x^5 \right ]_{-h}^h = 2 f_0 h + \frac{2}{3} b h^3 + \frac{2}{5}d h^5$

(3)

From (1),

$f_{-1} = f(x=-h) = f_0 - a h + b h^2 - c h^3 + d h^4$

(4)

and

$f_1 = f(x=h) = f_0 + a h + b h^2 + c h^3 + d h^4$

(5)

Therefore

$f_1+f_{-1} = 2 f_0 + 2 b h^2 + 2 d h^4$

(6)

$\frac{2}{3}bh^2 = \frac{1}{3}(f_1 + f_{-1}) - \frac{2}{3}f_0 - \frac{2}{3}d h^4$

(7)

Substituting into (2),

$I=\frac{h}{3}\left [ 4 f_0 + (f_1 + f_{-1}) - \frac{4}{5}d h^5 \right ]$

(8)

If we ignore the last term, this would create a truncation error proportional to the fourth derivative of the function f(x), which depends on $\inline h^5$ .

Consider a four times differentiable function $\inline f:I \subset \mathbb{R} \rightarrow \mathbb{R}$ , two distinct abscissas $\inline a < b \in I$ and a positive integer $\inline n$ . Then the following approximation holds:

$\int_a^b f(x) \mathrm{d}x \approx \frac{h}{3} \left[f(a) + f(b) + 4\sum_{k=1}^n f\left(a + (2k - 1)h\right) + 2\sum_{k=1}^{n-1} f(a + 2kh)\right],$

(9)

where $\inline \displaystyle h = \frac{b-a}{2n}$ , with an error bound of:

$\epsilon \leq \frac{(b-a)^5}{2880n^4} \sup_{x \in [a,b]} \left|\frac{d^4}{dx^4}f(x)\right|.$

(10)

This is a quadrature formula known as Simpson's rule. From the formula of the error bound you may notice that e.g. if one doubles the number of points n, the approximation error is decreased sixteen times. This rule approximates the area under f between a and b by the area under a quadratic polynomial P as seen in the following image, where m is the midpoint between a and b.

MISSING IMAGE!

1/SimpsonsRule-378.gif cannot be found in /users/1/SimpsonsRule-378.gif. Please contact the submission author.

Below you will find examples that show how to use each version of the function and the error estimates in each case. Generally this module provides better estimates than using the trapezoidal rule as implemented in Quadrature/trapezoidal. To convince yourself, compare the approximations generated by this module with the estimates given by the trapezoidal rule module.

References:

Mihai Postolache - "Metode Numerice", Editura Sirius
Wikipedia, http://en.wikipedia.org/wiki/Simpson's_rule

Example 1

References

Authors

Lucian Bentea (September 2006)

Simpson

doublesimpson(	int	`n`
	double	(`f`)(double)[function pointer]*
	double	`a`
	double	`b`	)

This version is to be used when you wish to pass the function to integrate as an argument, instead of passing an array with its values at distinct points.

Example 2

In what follows an approximation is found for the definite integral

$\int_1^3 \sin x \,dx$

(11)

and the absolute error from its actual value is estimated.

#include <codecogs/maths/calculus/quadrature/simpson.h>
#include <stdio.h>
#include <math.h>
 
// number of points
#define N 100
 
// function to integrate
double f(double x)
{
  return sin(x);
}
 
// the primitive of f, to estimate errors
double pf(double x)
{
  return -cos(x);
}
 
int main()
{
  // compute the approximate area
  double fi = Maths::Calculus::Quadrature::simpson(N, f, 1, 3),
 
  // use the Leibniz-Newton formula to find a more precise estimate
  realfi = pf(3) - pf(1); // uses the Leibniz-Newton formula
 
  // display problem data
  printf("      f(x) = sin(x)\n");
  printf("    points = %d\n\n", N);
 
  // display the result and error estimate
  printf("   I(1, 3) = %.15lf\n", fi);
  printf("real value = %.15lf\n", realfi);
  printf("     error = %.15lf\n\n", fabs(fi - realfi));
 
  return 0;
}

Output

f(x) = sin(x)
    points = 100
 
   I(1, 3) = 1.530294803828911
real value = 1.530294802468585
     error = 0.000000001360326

Parameters

n	the number of sample points of the function f, from which to approximate (must be even)
f	the function to integrate
a	the inferior limit of integration
b	the superior limit of integration

Returns

The definite integral of the given function from a to b.

Source Code

Source code is available when you buy a Commercial licence.

Not a member, then Register with CodeCogs. Already a Member, then Login.

Simpson

doublesimpson(	int	`n`
	const double*	`values`
	double	`a`
	double	`b`	)

This version is to be used when the function is not known analytically, but a table with its values at equally spaced abscissas is available.

Example 3

The code below gives an approximation for the definite integral

$\int_1^{1.3} \sqrt{x} \,dx$

(12)

and estimates the absolute error from its actual value.

#include <codecogs/maths/calculus/quadrature/simpson.h>
#include <stdio.h>
#include <math.h>
 
// the primitive of g, to estimate errors
double pf(double x)
{
  return 2*sqrt(x*x*x)/3;
}
 
int main()
{
  // values of the function at equally spaced abscissas
  double ordinates[7] = {1, 1.0247, 1.04881, 
  1.07238, 1.09544, 1.11803, 1.14017},
 
  // compute the approximate area
  fi = Maths::Calculus::Quadrature::simpson(7, ordinates, 1, 1.3),
 
  // use the Leibniz-Newton formula to find a more precise estimate
  realfi = pf(1.3) - pf(1); // uses the Leibniz-Newton formula
 
  // display problem data
  printf("      f(x) = sqrt(x)\n");
  printf("    points = 7\n\n");
 
  // display the result and error estimate
  printf(" I(1, 1.3) = %.15lf\n", fi);
  printf("real value = %.15lf\n", realfi);
  printf("     error = %.15lf\n\n", fabs(fi - realfi));
 
  return 0;
}