12,559,389 members (46,577 online)
alternative version

40.5K views
47 bookmarked
Posted

# Numerical Integration with Simpson's Rule

, 26 Jun 2014 BSD
 Rate this:
How to implement and use Simpson's rule

## Introduction

Simpson’s rule is a simple and effective technique for numerically evaluating integrals. However, practical implementation requires more than is often presented in introductions to the method. This article will show how to implement the algorithm n C# and will discuss how to prepare integrals for evaluation.

## Background

Calculus books often give a quick introduction to the trapezoid rule and Simpson’s rule as ways to approximately calculate integrals that cannot be computed exactly by hand.

Given an interval [a, b] and an even number n, Simpson’s rule approximates the integral

by

where h = (b - a)/n and xi = a + ih.

While that is a correct description of the algorithm, it leaves two related questions unanswered:

1. How good is the approximation?
2. How should I choose n?

In practice you seldom know before you start how large n needs to be in order to meet your accuracy requirements. If you can estimate your error at runtime, you can use that estimate to determine n.

A common approach is to repeatedly apply Simpson’s rule, doubling the number of integration points at each stage. The difference between two successive stages gives an estimate of the error, so you continue doubling the stages until the error estimate is below your tolerance.

There are a couple practical matters to implementing this approach. First, if naively implemented, the method will re-evaluate the integrand at the same points multiple times. At each new stage, half of the integration points are the same as they were at the previous stage. A little algebra makes it possible to reuse the integral estimate from the previous stage and only evaluate the integrand at the new integration points.

The second practical implementation matter is to not check the stopping condition too soon. If you check the stopping rule after each stage, it is possible in the early stages that your integration points miss regions of variability in your integrand and you are led to believe the method has converged when it has not. The code given here only checks for convergence after five stages of Simpson’s rule.

## Using the Code

The class SimpsonIntegrator has only two methods, overloaded implementations of Integrate. The simpler version of Integrate returns nothing but the value of the integral. The more detailed version returns the number of function evaluations used and an estimate of the error. The simpler version is implemented by calling the more detailed method and discarding the extra information. Also, the more detailed version allows the user to specify the maximum number of integration stages; the simpler version sets a default value for this parameter.

Here is an example integrating the function f(x) = cos(x) + 1 over the interval [0, 2π]:

double twopi = 2.0*Math.PI;
double epsilon = 1e-8;
double integral = SimpsonIntegrator.Integrate(x => Math.Cos(x) + 1.0, 0,
twopi, epsilon);
double expected = twopi;
Assert.IsTrue(Math.Abs(integral - expected) < epsilon);


Note that C#'s support for implicit functions is handy for simple integrands.

According to the theoretical derivation of Simpson’s rule, the method should integrate polynomials up to degree three exactly with only one stage of integration. Here's an example of calling the more detailed interface to specify that the method should use only one stage to integrate f(x) = x3 over an interval [a, b]. We also verify that the code respected our request to only use one stage by checking that the integrand is only evaluated three times. (The first stage of Simpson’s rule evaluates the integrate at each end of the interval and at the midpoint.)

double a = 2, b = 5;
double epsilon = 1e-10;
double integral;
int functionEvalsUsed;
double estimatedError;
double expected;
int maxStages = 1;

integral = SimpsonIntegrator.Integrate(x => x*x*x, a, b, epsilon,
maxStages, out functionEvalsUsed, out estimatedError);
expected = (b*b*b*b - a*a*a*a)/4.0;
Assert.AreEqual(functionEvalsUsed, 3);
Assert.IsTrue(Math.Abs(integral - expected) < epsilon);


Simpson’s rule only applies to smooth functions on a closed interval. If this does not describe your integral, it may still be possible to apply the method after a little preparation. On the other hand, naively applying the method directly where it is not appropriate can give poor results.

A function with a kink in the middle is not smooth and so Simpson’s rule does not directly apply. However, such a function can simply divided into two integrals where the integrand is smooth over both separately. For example, to integrate the function exp(-|x|) over [-3, 5], we would need to split the problem into one integral over [-3, 0] and another over [0, 5] because the integrand has a sharp bend at 0.

An integral over an infinite range cannot be computed by Simpson’s rule directly. However, it is often possible to use a change of variables to convert the integral into a well-behaved integral over a finite range.

Simpson’s rule also does not apply to functions that become infinite. See Care and treatment of singularities for an example of how to use Simpson’s rule for such integrals.

## Share

 Singular Value Consulting United States
I am an independent consultant in software development and applied mathematics. I help companies learn from their data to make better decisions.

Check out my blog or send me a note.

## You may also be interested in...

 First Prev Next
 Dynamic function parsing for SimpsonIntegrator Chuck G25-Sep-14 12:25 Chuck G 25-Sep-14 12:25
 Multiple integrals Alexandr Stefek28-Jun-14 10:11 Alexandr Stefek 28-Jun-14 10:11
 My vote of 5 Luis Fernando Pinilla29-Jun-13 10:28 Luis Fernando Pinilla 29-Jun-13 10:28
 epsilon Member 92430416-Mar-13 11:07 Member 9243041 6-Mar-13 11:07
 My vote of 5 David Serrano Martínez15-Feb-13 10:02 David Serrano Martínez 15-Feb-13 10:02
 Re: My vote of 5 Member 92430416-Mar-13 22:33 Member 9243041 6-Mar-13 22:33
 Re: My vote of 5 David Serrano Martínez7-Mar-13 8:38 David Serrano Martínez 7-Mar-13 8:38
 Epsilon is an arbitrarily small positive number used to control the integration error. Simpson`s numerical integration is based on partitions. If you want to integrate f(x) from A to B, you first split the segment A->B in N identical parts. Next you apply the algorithm evaluating f(x) on all the points of that partition and you get a result RESULT_1 (for instance). But if you split A->B with (say) 2*N identical parts, you get a better approach RESULT_2 (with respect to the true, unknown solution). You could next apply a finer partition (4*N), but... when to stop? If epsilon is small enough and the absolute difference abs(RESULT_2 - RESULT_1) is smaller than epsilon, then it is not worth doing your partition finer and finer any longer; there will no be gain for RESULT_3, RESULT_4, etc. And because you know Simpson's rule converges to the true value, epsilon is a good estimate of the maximum possible error you are getting. The second part of your comment is much trickier. In order to allow users to enter their own functions, you need a mathematical expression parser. Many math-specialized applications utilize them, but they are quite complicated tools.
 Simple and Works Member 818270720-Dec-12 7:09 Member 8182707 20-Dec-12 7:09
 Dum Dum Here Patrick Harris22-Feb-12 17:23 Patrick Harris 22-Feb-12 17:23
 Last Visit: 31-Dec-99 18:00     Last Update: 27-Oct-16 18:41 Refresh 1