13,595,454 members
alternative version

#### Stats

49.2K views
49 bookmarked
Posted 22 Feb 2012
Licenced BSD

# Numerical Integration with Simpson's Rule

, 26 Jun 2014
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
 Thanks Cryptonite28-Mar-17 6:07 Cryptonite 28-Mar-17 6:07
 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
 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
 My subject may not be the best to get attention, but it probably fits. Can you give an example or illustration of how this rule can be applied? I'm curious if I might learn something from this technique you offer to accomplish re-writing a project I took on at my job where I calculate the number of wave patterns, spaced at a given length, that fit into the circumference of a circle. Picture a sort of oscilloscope pattern if you will. To accomplish this, in the time frame I was given, I had to fudge some numbers. Based on that knowledge, it's hard to believe I accurately accomplished the goal (in one night), but I did. Needless to say, I would like to re-visit that project so I can run my calculations more mathematically and there by be able to enhance it more efficiently some day.
 Last Visit: 31-Dec-99 18:00     Last Update: 22-Jun-18 4:16 Refresh 1