Simpson's rule

Simpson's rules are methods used in mathematics to estimate the area under a curve. They are named after Thomas Simpson, who lived from 1710 to 1761. These rules help us find the value of definite integrals, which are important in many areas of science and engineering.

The most common rule is called Simpson's 1/3 rule, or simply Simpson's rule. It provides a way to approximate the area under a curve between two points. This rule works very well for curves that are smooth and can be exactly calculated for polynomials up to the third degree.

There is also a Simpson's 3/8 rule, which needs one more point to estimate the area. This rule can give more accurate results for certain kinds of curves. Both of these rules are special cases of a larger group of formulas called Newton–Cotes formulas, which are used in numerical analysis to solve complex mathematical problems.

Simpson's 1/3 rule

Simpson's 1/3 rule, also simply called Simpson's rule, is a method for numerical integration proposed by Thomas Simpson. It is based upon a quadratic interpolation and is the composite Simpson's 1/3 rule evaluated for n = 2. Simpson's 1/3 rule is as follows:

Since the error term is proportional to the fourth derivative of f at ξ, this shows that Simpson's rule provides exact results for any polynomial f of degree three or less, since the fourth derivative of such a polynomial is zero at all points.

If the second derivative f ″ exists and is convex in the interval ( a , b ), then

Derivations

Quadratic interpolation

Consider finding the area under a general parabola y = a x² + b x + c between x = − h and x = h for some positive number h. The midpoint of this interval is therefore at x = 0.

The area under the parabola, A, is therefore

A = ∫ − h h a x² + b x + c dx = [ a x³ / 3 + b x² / 2 + c x ] x = − h x = h = ( a h³ / 3 + b h² / 2 + c h ) − ( − a h³ / 3 + b h² / 2 − c h ) = 2 a h³ / 3 + 2 c h = h / 3 ( 2 a h² + 6 c ) .

Assuming the parabola has midpoint ( 0 , y₁ ) and endpoints ( − h , y₀ ) and ( h , y₂ ), substituting these three points into the parabola formula gives

y₀ = a h² − b h + c

y₁ = c

y₂ = a h² + b h + c .

Solving these gives

c = y₁

from the second equation, and

2 a h² = y₀ − 2 y₁ + y₂

by adding the first and third equations. Substituting this into the expression for A gives

A = h / 3 ( 2 a h² + 6 c ) = h / 3 ( y₀ − 2 y₁ + y₂ + 6 y₁ ) = h / 3 ( y₀ + 4 y₁ + y₂ ) .

Simpsons 1/3 rule approximates a definite integral over an interval [a, b] by replacing the integrand with a parabola which interpolates the function at x = a , x = b and the midpoint of the interval, which gives

∫ a b f ( x ) dx ≈ A = 1 / 3 h [ f ( a ) + 4 f ( a + h ) + f ( b ) ] .

Because of the 1 / 3 factor, Simpson's rule is also referred to as "Simpson's 1/3 rule" (see below for generalization).

Averaging the midpoint and the trapezoidal rules

Another derivation constructs Simpson's rule from two simpler approximations. For functions that behave like polynomials over the interval [ a , b ], the midpoint rule says M = ( b − a ) f ( a + b / 2 ) = ∫ a b f ( x ) dx − 1 / 24 ( b − a )³ f ″ ( a ) + O ( ( b − a )⁴ ) and the trapezoidal rule says T = ( b − a ) ( f ( a ) + f ( b ) / 2 ) = ∫ a b f ( x ) dx + 1 / 12 ( b − a )³ f ″ ( a ) + O ( ( b − a )⁴ ), where O ( ( b − a )⁴ ) denotes a term asymptotically proportional to ( b − a )⁴ . The two O ( ( b − a )⁴ ) terms are not equal; see Big O notation for more details. It follows from the above formulas that the leading error term vanishes if we take the weighted average 2 M + T / 3 = ∫ a b f ( x ) dx + O ( ( b − a )⁴ ). This weighted average is exactly Simpson's rule.

Using another approximation (for example, the trapezoidal rule with twice as many points), it is possible to take a suitable weighted average and eliminate another error term. This is Romberg's method.

Undetermined coefficients

The third derivation starts from the ansatz 1 b − a ∫ a b f ( x ) dx ≈ α f ( a ) + β f ( a + b / 2 ) + γ f ( b ) .

The coefficients α, β and γ can be fixed by requiring that this approximation be exact for all quadratic polynomials. This yields Simpson's rule. (This derivation is essentially a less rigorous version of the quadratic interpolation derivation, where one saves significant calculation effort by guessing the correct functional form.)

Composite Simpson's 1/3 rule

If the interval of integration [ a , b ] is in some sense "small", then Simpson's rule with n = 2 subintervals will provide an adequate approximation to the exact integral. By "small" we mean that the function being integrated is relatively smooth over the interval [ a , b ]. For such a function, a smooth quadratic interpolant like the one used in Simpson's rule will give good results.

However, it is often the case that the function we are trying to integrate is not smooth over the interval. Typically, this means that either the function is highly oscillatory or lacks derivatives at certain points. In these cases, Simpson's rule may give very poor results. One common way of handling this problem is by breaking up the interval [ a , b ] into n > 2 small subintervals. Simpson's rule is then applied to each subinterval, with the results being summed to produce an approximation for the integral over the entire interval. This sort of approach is termed the composite Simpson's 1/3 rule, or just composite Simpson's rule.

Suppose that the interval [ a , b ] is split up into n subintervals, with n an even number. Then, the composite Simpson's rule is given by

Dividing the interval [ a , b ] into n subintervals of length h = ( b − a ) / n and introducing the points x_i = a + i h for 0 ≤ i ≤ n (in particular, x₀ = a and x_n = b ), we have ∫ a b f ( x ) dx ≈ 1 / 3 h ∑ i = 1 n / 2 [ f ( x 2 i − 2 ) + 4 f ( x 2 i − 1 ) + f ( x 2 i ) ] = 1 / 3 h [ f ( x 0 ) + 4 f ( x 1 ) + 2 f ( x 2 ) + 4 f ( x 3 ) + 2 f ( x 4 ) + ⋯ + 2 f ( x n − 2 ) + 4 f ( x n − 1 ) + f ( x n ) ] = 1 / 3 h [ f ( x 0 ) + 4 ∑ i = 1 n / 2 f ( x 2 i − 1 ) + 2 ∑ i = 1 n / 2 − 1 f ( x 2 i ) + f ( x n ) ] .

This composite rule with n = 2 corresponds with the regular Simpson's rule of the preceding section.

The error committed by the composite Simpson's rule is − 1 / 180 h⁴ ( b − a ) f ( 4 ) ( ξ ) , where ξ is some number between a and b, and h = ( b − a ) / n is the "step length". The error is bounded (in absolute value) by 1 / 180 h⁴ ( b − a ) max ξ ∈ [ a , b ] | f ( 4 ) ( ξ ) | .

This formulation splits the interval [ a , b ] in subintervals of equal length. In practice, it is often advantageous to use subintervals of different lengths and concentrate the efforts on the places where the integrand is less well-behaved. This leads to the adaptive Simpson's method.

Examples

Approximating the natural logarithm of 2

Since

∫ 1 2 1 x dx = ln ( 2 ) − ln ( 1 ) = ln ( 2 ) ,

approximations of ln ( 2 ) can be generated by approximating this integral. Applying Composite Simpson's 1/3 rule with n = 6 intervals gives ln ( 2 ) = ∫ 1 2 1 x dx ≈ 1 / 18 ( 1 + 24 / 7 + 3 / 2 + 8 / 3 + 6 / 5 + 24 / 11 + 1 / 2 ) ≈ 0.69316 ,

which has a relative error of about 0.003 % .

An application to statistics

In statistics, when data tends around a central value with no left or right bias, it's said to be normally distributed. In the case where the mean is zero and the standard deviation is 1, the curve is said to follow the standard normal (or standard Gaussian) distribution. The equation of this distribution is

f ( x ) = 1 / √ 2 π exp ( − x² / 2 ) .

By the 68–95–99.7 rule, approximately 68.27% of values are within a single standard deviation of the mean, so

∫ − 1 1 1 / √ 2 π exp ( − x² / 2 ) dx ≈ 0.6827.

This result can be verified with Composite Simpson's 1/3 rule - applying the rule with n = 6 intervals gives ∫ − 1 1 1 / √ 2 π exp ( − x² / 2 ) dx ≈ 1 / √ 2 π × 1 / 9 ( 1 / e^{1 / 2} + 4 / e^{2 / 9} + 2 / e^{1 / 18} + 4 / 1 + 2 / e^{1 / 18} + 4 / e^{2 / 9} + 1 / e^{1 / 2} ) ≈ 1 / √ 2 π × 1.71142 ≈ 0.6827 ,

as expected.

Similarly, the 68–95–99.7 rule says approximately 95.45% of values are within two standard deviations of the mean, so

∫ − 2 2 1 / √ 2 π exp ( − x² / 2 ) dx ≈ 0.9545.

As before, this result can be verified with Composite Simpson's 1/3 rule - applying the rule with n = 6 intervals gives ∫ − 2 2 1 / √ 2 π exp ( − x² / 2 ) dx ≈ 1 / √ 2 π × 2 / 9 ( 1 / e² + 4 / e^{8 / 9} + 2 / e^{2 / 9} + 4 / 1 + 2 / e^{2 / 9} + 4 / e^{8 / 9} + 1 / e² ) ≈ 1 / √ 2 π × 2.39167 ≈ 0.9541 ,

which has a relative error of about 0.0419 % . The interval of integration can be shortened (thereby reducing discretization error) by noting that the integrand is an even function, so

∫ − 2 2 1 / √ 2 π exp ( − x² / 2 ) dx = 2 ∫ 0 2 1 / √ 2 π exp ( − x² / 2 ) dx .

Once again, applying Composite Simpson's 1/3 rule with n = 6 intervals gives 2 ∫ 0 2 1 / √ 2 π exp ( − x² / 2 ) dx ≈ 2 / √ 2 π × 1 / 9 ( 1 + 4 / e^{1 / 18} + 2 / e^{2 / 9} + 4 / e^{1 / 2} + 2 / e^{8 / 9} + 4 / e^{25 / 18} + 1 / e² ) ≈ 2 / √ 2 π × 1.19626 ≈ 0.9544 ,

which has an improved relative error of about 0.0104 % but with the same number of intervals.

Approximating π

Since

∫ 0 1 1 / 1 + x² dx = arctan ( 1 ) − arctan ( 0 ) = π / 4 ,

this can be rearranged to give

π = 4 ∫ 0 1 1 / 1 + x² dx .

Therefore, approximations of π can be generated by approximating this integral. Applying Composite Simpson's 1/3 rule with n = 6 intervals gives π = 4 ∫ 0 1 1 / 1 + x² dx ≈ 2 / 9 ( 1 + 144 / 37 + 9 / 5 + 16 / 5 + 18 / 13 + 144 / 61 + 1 / 2 ) ≈ 3.141591 ,

which remarkably only has a relative error of about 0.00002 % .

Determining the number of intervals for a desired accuracy

Suppose we wish to determine the number of intervals required to approximate ∫ 0 π sin ( x ) dx with an absolute error of less than 0.00001 . The error term in Composite Simpson's 1/3 rule is

− π h⁴ / 180 sin ( ξ )

for some ξ between 0 and π . Since the absolute error is to be less than 0.00001 , we can calculate

so n = 22 will generate the required accuracy.

For comparison purposes, suppose we wish to be assured of this degree of accuracy using the Composite Trapezoidal Rule. In this case, the error term is

− π h² / 12 sin ( ξ )

for some ξ between 0 and π . Since the absolute error is to be less than 0.00001 , we can calculate

so n = 509 will guarantee the required accuracy. This is significantly more calculations compared to Composite Simpson's 1/3 rule.

Simpson's 3/8 rule

Simpson's 3/8 rule is another way to estimate areas under curves, named after Thomas Simpson. It uses a special method that fits a curved shape to four points on the curve, making it more accurate for certain problems.

This rule works by dividing the area into parts and using a formula that involves adding up values from these points. It’s especially useful when you need better accuracy than simpler methods, though it needs one more point to do so. This makes it a handy tool for solving complex math problems that can’t be solved exactly.

Alternative extended Simpson's rule

This is another way to use Simpson's rule for estimating areas under curves. Instead of looking at separate pieces, this method looks at pieces that overlap. The formula uses different weights for different points along the curve to get a better estimate.

Simpson's rules can also help estimate the area under very narrow peak shapes. These special rules use points both inside and just outside the area being measured to improve accuracy. These methods are linked to a mathematical idea called the Euler–MacLaurin formula, which helps make these estimates even better.

Composite Simpson's rule for irregularly spaced data

When we want to find the area under a curve, sometimes the data points we have are not evenly spaced. This can happen if we missed some points or if they were unevenly collected.

To handle this, we can use a version of Simpson's rule that works with uneven intervals. If we split the area into an even number of parts, we can use a special formula to estimate the area. If there are an odd number of parts, we use the formula for the first part and then add a bit more for the last part. This helps us get a good estimate even when the data points are not evenly spread out.

Numerical stability

Simpson's Rule is known for being stable with small mistakes that can happen when doing calculations. This stability is shared by all Newton–Cotes formulas.

When using Simpson's Rule on a function over a range, small mistakes in each step add up. However, these total mistakes stay small no matter how tiny the steps become. This makes Simpson's Rule much better than methods used for finding slopes, which can become very wrong with small changes.