Now Loading...

Please download on PC (not smartphone) and extract the ZIP file. Then double-click VCSSL.bat (batch file) to execute for Microsoft® Windows®, or execute VCSSL.jar (JAR file) on the command line for Linux® and other OSes.
For details, see How to Use.

Numerical Integration Using Simpson's Rule (Quadratic Approximation)

This VCSSL program is the third part of our series on numerical integration, continuing from the previous article, and again numerically computes the value of a definite integral.

This time, we go a step beyond the previous methods by using Simpson's Rule, which achieves even higher accuracy when summing up the areas of small regions under the curve. Although the difference in code from the rectangular method and the trapezoidal method comes down to just a single line, the accuracy improvement is dramatic.

Please note that this program is a sample implementation for explaining the algorithm, so it does not come with a fancy GUI for parameter input. To modify the function being integrated or change other parameters, please edit the relevant parts of the source code directly.

The code is written in VCSSL, but feel free to translate it into other languages as needed (a C version and a C++ version are also included below). The code is released under a public domain / CC0 license, so you are free to use or modify it at your own discretion.

[ Articles in the Numerical Integration Series ]

• Part 1: Numerical Integration Using the Rectangular Method
• Part 2: Numerical Integration Using the Trapezoidal Method (previous)
• Part 3: Numerical Integration Using Simpson's Rule (this article)
• Summary: General Numerical Integration Programs

How to Use

Download and Extract

At first, click the "Download" button at the above of the title of this page by your PC (not smartphone). A ZIP file will be downloaded.

Then, please extract the ZIP file. On general environment (Windows®, major Linux distributions, etc.), you can extract the ZIP file by selecting "Extract All" and so on from right-clicking menu.

» If the extraction of the downloaded ZIP file is stopped with security warning messages...

Execute this Program

Next, open the extracted folder and execute this VCSSL program.

For Windows

Double-click the following batch file to execute:

For Linux, etc.

Execute "VCSSL.jar" on the command-line terminal as follows:

» If the error message about non-availability of "java" command is output...

After Launching

Once the program launches, it runs the numerical integration process and displays a graph of the integrated function. The console also outputs the calculated integral values.

Modifying the Function and Parameters

If you want to change the function, integration range, or precision, open the program file named IntegralSimpson.vcssl in a text editor and modify the relevant parts of the code.

Topic Overview

Recap: Rectangular and Trapezoidal Methods

This article builds on the previous installments, so let's begin with a quick recap. If you'd like to review the foundation more thoroughly, feel free to revisit the previous articles:

• Part 1: Numerical Integration Using the Rectangular Method
• Part 2: Numerical Integration Using the Trapezoidal Method

In Part 1, we computed the integral by approximating the area under the curve of \(y = f(x)\) using thin rectangular strips. Here's a visual representation:

Concept of the Rectangular Method -- Approximating the area under \(y = f(x)\) using narrow rectangles

Since the area under the curve \(f(x)\) represents the value of the integral, summing up the areas of these rectangles gives an approximate result.

As you may recall, the area between the function \(f(x)\) and the x-axis theoretically equals the integral of \(f(x)\). So by computing that area numerically in a program, you can approximate the integral -- even without solving the integral symbolically.

In the rectangular method, we simply summed the areas of those thin rectangles. But this approach introduced noticeable error, which is easy to understand from the image above: The upper edge of the region should be a smooth curve, but the rectangles create a jagged approximation, leading to either overestimation or underestimation of the actual area. That was the source of the error.

To address this, in Part 2 we switched to the trapezoidal method, replacing each rectangle with a **trapezoid**, which better conforms to the curve and reduces error. Here's what that looks like:

Concept of the Trapezoidal Method -- Approximating the area with narrow trapezoids

By summing the areas of trapezoids instead of rectangles, we can achieve a more accurate estimate of the integral.

Just by looking at the image, you can see that the jagged top edges of the rectangles are smoothed out by the trapezoids -- intuitively suggesting less error.

And in fact, when we decreased the strip width \(\Delta x\), the rectangular method showed a linear improvement in accuracy, while the trapezoidal method showed a quadratic improvement.

For example, reducing \(\Delta x\) by a factor of 10 cut the error by a factor of 10 with rectangles, but by a factor of 100 with trapezoids. With \(\Delta x\) reduced to one-hundredth, trapezoids gave us one ten-thousandth the error that rectangles did. That's a big difference.

How Can We Improve Even Further? -- Approximate the Top Edge with a Curve

Now that we've reviewed the basics, let's move on to the main topic. We'll make another small improvement to go beyond the trapezoidal method.

Don't worry -- it's not a major overhaul. Just like before, all we're doing is changing one line in the code that computes the area of a strip. So once again, everything boils down to: "How precisely can we estimate the area of each tiny strip?"

Let's take a closer look at the source of the error in the trapezoidal method. Here's one of those narrow strips (with the curvature exaggerated for clarity):

A single narrow strip, zoomed in

When using trapezoids, the top edge approximates the curve \(f(x)\) with a straight line. This introduces a small error -- the orange regions -- because the area of the trapezoid doesn't exactly match the true area under the curve.

As you can see, although it's better than the rectangles, the top edge is still a straight line. This means we're approximating the curved function \(f(x)\) with a straight segment, causing a mismatch -- the orange areas in the diagram.

Depending on the curvature, this leads to small over- or under-estimations, which add up across the integration range and become the total error.

So what can we do? The key idea is simple: If the top edge of the strip is curved, why not approximate it with a curve instead of a straight line?

Using a curve to approximate the top edge

If we replace the straight line with a curve, it intuitively feels like the error will decrease, and the result will be more accurate.

Just by looking at this, it seems clear that we'll reduce the over- and under-estimations, leading to better accuracy.

What Kind of Curve Should We Use to Approximate the Top of Each Strip? -- Use a Simple Quadratic Function That Can Be Integrated by Hand

Now we reach the key question: what kind of curve should we use to approximate \(f(x)\) along the top edge of each small strip? This choice is at the heart of Simpson's Rule.

In theory, the most accurate approach would be to just use the original function \(f(x)\) itself. But if we did that, we'd need to integrate \(f(x)\) by hand to get the area of each strip -- which defeats the purpose. After all, we're writing a program to handle integration precisely because integrating \(f(x)\) by hand may be difficult or even impossible.

So what do we use instead? A quadratic function. Quadratic curves are simple, as shown in the graph below, and as you may remember from high school math, they can be integrated exactly by hand.

So, if we can approximate \(f(x)\) with a well-chosen quadratic, then the area of each small strip can be computed with a clean formula -- and we can just replace the trapezoid-area formula from the previous method with this new one.

Various Graphs of Quadratic Functions

The shape and position of a quadratic curve change depending on its coefficients.

As you might recall from school, the shape and position of a quadratic function depend on its coefficients (see the figure above). Maybe you've wrestled with graphing something like \(y = 2 x^2 + 3 x - 4\) for a test.

Here, our goal is to approximate \(f(x)\), so we'll choose the coefficients of our quadratic so that it hugs the shape of \(f(x)\) as closely as possible within each small strip.

There are many ways to choose the coefficients of a polynomial to fit a function. But wefd like to keep things simple.

So here we'll go with a straightforward approach: we'll choose the coefficients so that the quadratic passes through three points: the left end, the right end, and the midpoint of the strip.

If you only fit the endpoints, there are infinitely many possible parabolas, some curving upward and some downward. By requiring it to also pass through the midpoint, we can ensure the curve better fits the actual shape of \(f(x)\).

Approximating the Top of Each Strip with a Quadratic Curve

The coefficients of the quadratic function are chosen so that it passes through the same three points as \(f(x)\): the left end, right end, and midpoint of the strip. The resulting area can be expressed using Simpson's Rule.

Once we choose the coefficients so that the quadratic passes through the three specified points, we've locked in one unique curve. (We'll skip the manual derivation here -- it can be done with pencil and paper if you like.)

Then, we integrate that quadratic to get the area of the strip. This yields a famous formula known as Simpson's Rule, and it looks like this:
\[ \textrm{Area of i-th strip} = \frac{ f(x_i) + f(x_i + \Delta x) + 4f(x_i+ \Delta x / 2) }{6} \Delta x \]
Here, \( x_i \) is the left endpoint of the strip, \( x_i+\Delta x \) is the right endpoint, \( x_i + \Delta x/2 \) is the midpoint, and \( \Delta x \) is the strip width.

Now just sum this over all the strips:
\[ \textrm{Integral} = \sum_{i=0}^{N-1} \frac{ f(x_i) + f(x_i + \Delta x) + 4f(x_i+ \Delta x / 2) }{6} \Delta x \]

You can rearrange this a bit to avoid duplicate evaluations, but this version maps most clearly to the code, so we'll keep it in this form.

Basically, we've just swapped out the area formulas we used in the rectangular and trapezoidal methods with a more accurate one based on a quadratic approximation.

Code

Now let's take a look at how to implement numerical integration using Simpson's Rule in actual code. The following program is written in VCSSL, which has a syntax that should feel quite familiar if you're used to C-like languages.

For reference, C and C++ versions of the code are also included below.

Full Code

Here's the full VCSSL implementation. It's a short program with fewer than 50 lines of code:

coding UTF-8;         // Specify character encoding (to avoid text corruption)
import Math;          // Import math function library
import tool.Graph2D;  // Import graphing library


// Parameter constants for integration
const double A = 0.0;  // Lower bound of integration interval
const double B = 1.0;  // Upper bound of integration interval
const int N = 10;      // Number of small intervals (larger = finer = more accurate)

// The function to be integrated: f(x) = cos(x)
double f(double x) {
	return cos(x);
}

// ---------------------------------------------
// Numerical integration using Simpson's Rule
// ---------------------------------------------

double delta = (B - A) / N;  // Width of a small interval (Delta x)
double value = 0.0;          // Result of integration (total area)

for(int i=0; i<N; i++) {
	
	// X-coordinate of the left edge of the i-th interval
	double x = A + i * delta;
	
	// Output the current integration result up to x
	println(x, value); // (Comment out this line for large N to improve speed)
	
	// Add the area of the i-th interval to the total
	value += ( f(x) + f(x+delta) + 4.0 * f(x+delta/2.0) ) * delta / 6.0;
}

// Output the final result of the integration
println(B, value);

// Plot the output data to a graph
// (Comment out this line for large N to avoid slowdown)
setGraph2DData(newGraph2D(), load());

That's all. The general flow is:

• Define the target function and parameters at the beginning,
• Use a "for" loop in the middle to perform numerical integration,
• And finally, visualize the integration process as a graph.

The C version of the code looks like this:

The C++ version of the code is shown below:

That's it for the C and C++ versions. Next, let's take a closer look at how it works in detail.

Just One Line Changed from the Previous Version!

The details of the code structure haven't changed much from the previous articles on the rectangular and trapezoidal methods, so for those explanations, please refer to the code walkthrough for the rectangular method.

Compared to the code for the rectangular and trapezoidal methods, there's only one line that's different here.

[ ↓ Rectangular Method Code ]

	...
	
	// Add the area of the i-th strip (rectangle)
	value += f(x) * delta;
	
	...

[ ↓ Trapezoidal Method Code ]

	...
	
	// Add the area of the i-th trapezoid
	value += ( f(x) + f(x+delta) ) * delta / 2.0;
	
	...

[ ↓ This Time (Simpson's Rule) ]

	...
	
	// Add the area of the i-th small region
	value += ( f(x) + f(x+delta) + 4.0 * f(x+delta/2.0) ) * delta / 6.0;
	
	...

As you can see, we're simply replacing the area calculation using rectangles or trapezoids with the formula based on Simpson's Rule introduced earlier, which uses a quadratic approximation:
\[ \textrm{Area of i-th strip} = \frac{ f(x_i) + f(x_i + \Delta x) + 4f(x_i+ \Delta x / 2) }{6} \Delta x \]

So, How Much Error Is Reduced?

Now let's see how much the error is reduced by this single-line change. When you run this program, the integration process is visualized on a graph.

Since the target function is currently set to \(f(x) = \cos x\), the output graph should resemble the shape of \(\sin x\), which is the integral of \(\cos x\).

A graph plotting the integration result over time

The red line (almost invisible because it overlaps the blue one) shows the integral values printed by the program and plotted directly. The blue line is a plot of the exact integral function \(sin x\).

As shown in the figure above, the red line (hard to see due to overlap) is the output of our code using Simpson's Rule, while the blue line represents the exact value of the integral, \(sin x\), for comparison.

Even with the trapezoidal method in the previous article, they overlapped quite closely -- but this time, they're so close it's almost impossible to distinguish them. This is despite the fact that we're using only \(N = 10\), which is a relatively coarse subdivision.

If you zoom in really far, you might see a slight gap between the calculated results and the true curve. But at this level, it's easier to examine the difference by looking directly at the numerical output in the console, which might look like this:

In this table, the left column shows the x-values as we integrate from a to b, and the right column shows the cumulative integral values at each point. (This is the same data used to generate the graph above.)

Each row adds the area of one small region to the total integral. The final value, 0.8414710140343371, is the result of:
\[ \int^{1}_{0} \cos(x) dx \]

The exact answer is \(\sin 1 \approx 0.8414709848078965...\) so they differ starting from the 7th digit after the decimal -- that's our error.

With the trapezoidal method in the previous article, the difference appeared at the 4th digit. With the rectangular method, the difference started at the 3rd digit.

As we increase the number of divisions \(N\) (which makes the step width \(\Delta x = (b-a)/N\) smaller), the error should become smaller -- this was discussed in more detail in the rectangular method article.

Let's try actually increasing \(N\) in the program and check the final integral value:

To make it easier to see, we've highlighted the digits that differ from the exact value \(\sin 1\).

From \(N = 1\) to \(N = 1000\), each 10× increase in \(N\) boosts the number of matching digits by roughly 4 digits -- aside from the last 1-2 digits, which behave slightly differently (more on that below). This confirms that Simpson's Rule is a 4th-order method.

For comparison:

• The trapezoidal method is 2nd-order, improving by about 2 digits per 10× increase in \(N\).
• The rectangular method is 1st-order, improving by only 1 digit per 10× increase in \(N\).

So with Simpson's Rule, a small change in \(N\) can lead to a huge improvement in precision -- and all this comes from a single-line difference in the code.

This kind of dramatic leap is one of the most exciting things about computational algorithms.

Even when a method seems complex at first glance, actually writing the code yourself and seeing it output highly accurate results is a thrilling and rewarding experience.

[ Articles in the Numerical Integration Series ]

• Part 1: Numerical Integration Using the Rectangular Method
• Part 2: Numerical Integration Using the Trapezoidal Method (previous)
• Part 3: Numerical Integration Using Simpson's Rule (this article)
• Summary: General Numerical Integration Programs

License

This VCSSL/Vnano code (files with the ".vcssl" or ".vnano" extensions) is released under the CC0 license, effectively placing it in the public domain. If any sample code in C, C++, or Java is included in this article, it is also released under the same terms. You are free to use, modify, or repurpose it as you wish.

* The distribution folder also includes the VCSSL runtime environment, so you can run the program immediately after downloading. The license for the runtime is included in the gLicenseh folder.
(In short, it can be used freely for both commercial and non-commercial purposes, but the developers take no responsibility for any consequences arising from its use.) For details on the files and licenses included in the distribution folder, please refer to "ReadMe.txt".

* The Vnano runtime environment is also available as open-source, so you can embed it in other software if needed. For more information, see here.