Numerical Integration

MAT 229, Spring 2025

Supporting materials

Stewart’s Calculus (Long’s overview)
	Section 7.7: Approximate integration

Strang’s Calculus:
	Vol 2, Sec 3.6: Numerical Integration

Boelkins/Austin/Schlicker’s Active Calculus:
	Section 5.6: Numerical integration

Overview

In this lab we’ll see a few different methods for computing the numerical approximations we’ve learned this week. 

One thing we’ll see is that we have another super-power: to create functions, even functions of several variables. 

Our key new commands this week:
	General commands:
		Clear  -- clears the definition of a variable, function, etc.
		Divide Cell
	Specific commands for the integration:
		Integrate, NIntegrate (to check our work)
		Rest -- all but the first
		Most -- all but the last
		Total -- add up everything
		Table -- build a list
		Sum -- summation
		Module	
	Options:	
		Plot (Option: AxesOrigin)

Rectangles and Partitions

Suppose that we want to approximate the area between the x-axis and y=(e^(-(x^2))) for 0≤x≤2 using four rectangles:

We need to break the interval up into a partition --  we need the x values, and then the corresponding y values. We usually use equally spaced points, and we will for our work here.

Suppose we are integrating from a to b, and using n rectangles: let

Δ x=((b-a)/n) and (Subscript[x, k])=a+k Δ x.

Thus the x values will run from (Subscript[x, 0]) to (Subscript[x, n]). We construct a list of those values.

One of the cool things about Mathematica (and many macro languages, frankly) is their ability to work with lists. For example, if we have a list of x values, call it xs, then we can  create the list of y-values ys by simply computing f[xs]. So f can apply to a number, or to a list. Very handy and  flexible!

We will start “the hard way”, and then build some better tools to quickly integrate any function using our various methods. I illustrated this in class, but now it’s your turn.... For this crude approach, we’ll use the commands Total, Table, Rest, and Most.

We’ll integrate function f over the interval [a,b], using 4 subintervals and all five of our methods.

f[x_]:=E^(-x^2)
Plot[f[x],{x,a,b},AxesOrigin->{0,0}]
a=0
b=2
integral=NIntegrate[f[x],{x,a,b}]

n=4
deltaX=(b-a)/n
xs=Table[xi,{xi,a,b,deltaX}]
ys=f[xs]

rrr=N[Total[Rest[ys]]*deltaX]
lrr=N[Total[Most[ys]]*deltaX]
trap=(lrr+rrr)/2

midpoints=Most[f[xs+0.5deltaX]]
mid=Total[midpoints]*deltaX
simp=(2mid+trap)/3
integral

More general methods for using these rules

Each of our basic rules (from which all five are constructed) can be written as a summation, as follows.

Left endpoint rule: (∫_a^b)f(x) d x≈(∑_(k=1)^n)f((Subscript[x, k-1]))Δ x=(f((Subscript[x, 0]))+f((Subscript[x, 1]))+…+f((Subscript[x, n-1])))Δ x
Right endpoint rule: (∫_a^b)f(x) d x≈(∑_(k=1)^n)f((Subscript[x, k]))Δ x=(f((Subscript[x, 1]))+f((Subscript[x, 2]))+…+f((Subscript[x, n])))Δ x
Midpoint rule: (∫_a^b)f(x) d x≈(∑_(k=1)^n)f((((Subscript[x, k-1])+(Subscript[x, k]))/2))Δ x=(f((((Subscript[x, 0])+(Subscript[x, 1]))/2))+f((((Subscript[x, 1])+(Subscript[x, 2]))/2))+…+f((((Subscript[x, n-1])+(Subscript[x, n]))/2)))Δ x

The trapezoidal rule for estimating (∫_a^b)f(x) d x with n trapezoids is the average of the left endpoint rule and the right endpoint rule each with n rectangles:

	(∫_a^b)f(x) d x≈((((f((Subscript[x, 0]))+f((Subscript[x, 1]))+f((Subscript[x, 2]))+…+f((Subscript[x, n-2]))+f((Subscript[x, n-1])))Δ x+(f((Subscript[x, 1]))+f((Subscript[x, 2]))+f((Subscript[x, 3]))+…+f((Subscript[x, n-1]))+f((Subscript[x, n])))Δ x)/2)
	=(1/2)Δ x(f((Subscript[x, 0]))+2f((Subscript[x, 1]))+…+2f((Subscript[x, n-1]))+f((Subscript[x, n])))
	
As you will recall, the trapezoidal rule is about half as accurate as the Midpoint method, in general.

Let’s build the left rectangle rule as an example

We will use the Module and Sum commands. We’re also going to build a local copy of the integrand function. 

The rule will be called like this: 
	leftRectangleRule[x^3,{x,0,1,2}]
or
	leftRectangleRule[f[x],{x,0,1,2}]
if you’ve defined f[x] (notice that I start my command name with a lower-case “l”). 

Once we build the left, you build rightRectangleRule! Then the trapezoidal will be easy: you can just use the average.

Clear[lrr,f,a,b,n,dx]
leftRectangleRule[f_,a_,b_,n_]:=
Module[
{dx=N[(b-a)/n]},
dx(Sum[f[a+k dx],{k,0,n-1}])
]
leftRectangleRule[x^3,0,1,2]
f[x_]:=x^3
leftRectangleRule[f,0,1,2]

leftRectangleRule[exp_,{x_,a_,b_,n_}]:=
Module[
{dx=N[(b-a)/n],f},
f[u_]:=exp/.x->u;
dx(Sum[f[a+k dx],{k,0,n-1}])
]
leftRectangleRule[x^3,{x,0,1,2}]
f[x_]:=x^3
leftRectangleRule[f[x],{x,0,1,2}]
leftRectangleRule[Sin[x],{x,0,1,2}]

Once we build the left, you build rightRectangleRule.

rightRectangleRule[exp_,{x_,a_,b_,n_}]:=
Module[
{dx=N[(b-a)/n],f},
f[u_]:=exp/.x->u;
dx(Sum[f[a+k dx],{k,1,n}])
]
rightRectangleRule[x^3,{x,0,1,2}]
f[x_]:=x^3
rightRectangleRule[f[x],{x,0,1,2}]
rightRectangleRule[Sin[x],{x,0,1,2}]

Then the trapezoidal will be easy: you can just use the average.

trapezoidalRule[exp_,{x_,a_,b_,n_}]:=
(leftRectangleRule[exp,{x,a,b,n}]+
rightRectangleRule[exp,{x,a,b,n}])/2
trapezoidalRule[x,{x,0,1,2}]
trapezoidalRule[x^3,{x,0,1,2}]
f[x_]:=x^3
trapezoidalRule[f[x],{x,0,1,2}]
trapezoidalRule[Sin[x],{x,0,1,2}]

Now we can reproduce the result of this example quite easily:

Estimate (∫_1^3)(1/(1+(x^2)))d x using the trapezoid rule with 4 trapezoids.

n=4, the number of trapezoids

Δ x=((3-1)/4)=(1/2)

Using (Subscript[x, k])=a+k Δ x with a=1 and Δ x=(1/2)=0.5

(Subscript[x, 0])=1 ((Subscript[x, 0]) is always the lower limit, the a-value.)

(Subscript[x, 1 ])=1.5

(Subscript[x, 2])=2

(Subscript[x, 3])=2.5

(Subscript[x, 4])=3 ((Subscript[x, n]) is always the upper limit, the b-value.)

(∫_1^3)(1/(1+(x^2)))d x≈(1/2)Δ x(f((Subscript[x, 0]))+2f((Subscript[x, 1]))+2f((Subscript[x, 2]))+2f((Subscript[x, 3]))+f((Subscript[x, 4])))
	=(1/2)((1/2))((1/(1+(((1))^2)))+2(1/(1+(((1.5))^2)))+2(1/(1+(((2))^2)))+2(1/(1+(((2.5))^2)))+(1/(1+(((3))^2))))≈0.472812

trapezoidalRule[1/(1+x^2),{x,1,3,4}]

trapezoid[exp_,{x_,a_,b_,n_}]:=Module[{dx=N[(b-a)/n],f},
f[u_]:=exp/.x->u;
Return[0.5dx(f[a]+f[b]+2Sum[f[a+k dx],{k,1,n-1}])]
]
midpoint[exp_,{x_,a_,b_,n_}]:=Module[{dx=N[(b-a)/n],f},
f[u_]:=exp/.x->u;
Return[dx Sum[f[a+k dx - 1/2dx],{k,1,n}]]
]
simpson[exp_,{x_,a_,b_,n_}]:=(trapezoid[exp,{x,a,b,n}]+2midpoint[exp,{x,a,b,n}])/3
trapezoid[1/(1+x^2),{x,1,3,4}]
midpoint[1/(1+x^2),{x,1,3,4}]
simpson[1/(1+x^2),{x,1,3,4}]
N[Integrate[1/(1+x^2),{x,1,3}]]

Questions

Use the fundamental theorem of calculus to get the exact value for (∫_1^3)(1/(1+(x^2)))d x.

Integrate[1/(1+x^2),{x,1,3}]
ArcTan[3]-ArcTan[1]

How good is the trapezoidal approximation we just got? Find the absolute error |exact-approximation|.

Abs[(ArcTan[3]-ArcTan[1])-0.4728116710875332]

Now let’s build the midpointRule.

midpointRule[exp_,{x_,a_,b_,n_}]:=
Module[
{dx=N[(b-a)/n],f,shift},
f[u_]:=exp/.x->u;
shift=a+dx/2;
dx(Sum[f[shift+k dx],{k,0,n-1}])
]
midpointRule[x^3,{x,0,1,2}]
f[x_]:=x^3
midpointRule[f[x],{x,0,1,2}]
midpointRule[Sin[x],{x,0,1,2}]

Example

Compare the trapezoid and midpoint rules to estimate (∫_0^2)(Sqrt[1+(x^3)])d x with 5 trapezoids

Then compute their Simpson’s Rule average. The Simpson’s Rule average would be  considered an n=10 rule, since it uses an eleven point partition.

Errors in numerical integration

When you approximate a quantity, typically you need to have some idea of how good your approximation is.

Definition

The absolute error E in estimating value V with approximation A is
	E=|V-A|.

Error bounds

Suppose |f''(x)|≤K for a≤x≤b where K is a constant. The absolute error in approximating (∫_a^b)f(x)d x with

the trapezoid rule using n trapezoids is ≤(((K(b-a))^3)/(12(n^2))).

the midpoint rule using n rectangles is ≤(((K(b-a))^3)/(24(n^2))).

Notice that the error bound for Midpoint is half that of the Trapezoidal rule (meaning that Midpoint is generally about twice as good as Trapezoidal).

Example

We want to approximate (∫_0^1)(e^(-(x^2)))d x. Let f(x)=(e^(-(x^2))).

What is f''(x)?

f'(x)=-2x (e^(-(x^2)))

f''(x)=-2(e^(-(x^2)))+4(x^2)(e^(-(x^2)))=2(e^(-(x^2)))(2(x^2)-1)

We want a constant value for K, in |f''(x)|≤K for 0≤x≤1. Plot the graph of |f''(x)| for these values of x.

Plot[Abs[2E^(-x^2)(2x^2-1)],{x,0,1}]

From this plot, choose K to be at least the maximum height. Here that is 2.
	K=2

Estimate (∫_0^1)(e^(-(x^2)))d x with the trapezoid rule using 4 trapezoids.

n=4, the number of trapezoids

Δ x=((1-0)/4)=(1/4)=0.25

Using (Subscript[x, k])=a+k Δ x with a=0 and Δ x=0.25

(Subscript[x, 0])=0

(Subscript[x, 1 ])=0.25

(Subscript[x, 2])=0.5

(Subscript[x, 3])=0.75

(Subscript[x, 4])=1

(∫_0^1)(e^(-(x^2)))d x≈(1/2)Δ x(f((Subscript[x, 0]))+2f((Subscript[x, 1]))+2f((Subscript[x, 2]))+2f((Subscript[x, 3]))+f((Subscript[x, 4])))
	=(1/2)((1/4))((e^(-(((0))^2)))+2(e^(-(((0.25))^2)))+2(e^(-(((0.5))^2)))+2(e^(-(((0.75))^2)))+(e^(-(((1))^2))))≈0.742984

trapezoidalRule[(E^(-(x^2))),{x,0,1,4}]

The absolute error in our approximation is
	|exact value-approximation|=|exact value-0.742984|
	≤(((K(b-a))^3)/(12(n^2)))=((2(((1-0))^3))/(12((4^2))))=0.0104167
The approximation is good to at least one decimal place. There might be just a little bit of error in the second decimal place.

Question

How many trapezoids should we use to make sure the trapezoid approximation to (∫_0^1)(e^(-(x^2)))d x has error less than 0.0001?

Simpson’s rule

Estimation formula

To approximate (∫_a^b)f(x)d x, let n be an even number and Δ x=((b-a)/n), then Simpson’s rule is
	(∫_a^b)f(x)d x≈((Δ x)/3)(f((Subscript[x, 0]))+4f((Subscript[x, 1]))+2f((Subscript[x, 2]))+4f((Subscript[x, 3]))+2f((Subscript[x, 4]))+…+2f((Subscript[x, n-2]))+4f((Subscript[x, n-1]))+f((Subscript[x, n])))

Simpson’s rule is actually just an average of the Trapezoidal rule and the Midpoint rule. Simpson’s 2n rule is the average of the Trapezoidal and Midpoint rules with n rectangles:
						(Subscript[S, 2n])=(((Subscript[T, n])+2(Subscript[M, n]))/3)
Notice that we’re adding two copies of the Midpoint rule to one copy of the Trapezoidal rule -- so we have three estimates of the integral (so we divide by 3). Because the error of the Trapezoidal rule tends to be about twice the error of the Midpoint rule, and of opposite sign, the errors cancel in a beautiful way to create a cool new and improved rule!

Our Simpson’s rule (simpsonsRule) will just be the weighted average of the midpoint and trapezoidal rules.

simpsonsRule[exp_,{x_,a_,b_,n_}]:=
(2midpointRule[exp,{x,a,b,n}]+
trapezoidalRule[exp,{x,a,b,n}])/3
simpsonsRule[x,{x,0,1,2}]
simpsonsRule[x^3,{x,0,1,2}]
f[x_]:=x^3
simpsonsRule[f[x],{x,0,1,2}]
simpsonsRule[Sin[x],{x,0,1,2}]

Error bound

The absolute error in the above approximation is
	error≤(((K(b-a))^5)/(180(n^4)))
where K is a constant such that |(f^((4)))(x)|≤K for a≤x≤b.

Remember that Simpson’s Rule 2n is built off of trap and mid with n subdivisions. This is important given the form of the error terms.

Example

We want to approximate (∫_0^1)(e^(-(x^2)))d x. Let f(x)=(e^(-(x^2))).

f[x_]=(E^(-(x^2)))
f''''[x]

The fourth derivative is
	(f^((4)))(x)=12 (e^(-(x^2)))-48 (e^(-(x^2))) (x^2)+16 (e^(-(x^2))) (x^4)
and the graph of its absolute value is

Plot[Abs[f''''[x]],{x,0,1}]

We want a constant value for K, where |(f^((4)))(x)|≤K. The maximum y-value from the above graph is 12. Use K=12 (or anything bigger than 12 -- but 12 gives the “tightest bound” on the error).

To estimate (∫_0^1)(e^(-(x^2)))d x with the Simpson’s rule using n=4:

Δ x=((b-a)/n)=((1-0)/4)=0.25

(Subscript[x, 0])=0
(Subscript[X, 1])=0.25
(Subscript[x, 2])=0.5
(Subscript[x, 3])=0.75
(Subscript[x, 4])=1

Using n=4, Simpson’s rule is
	(∫_a^b)f(x)d x≈((Δ x)/3)(f((Subscript[x, 0]))+4f((Subscript[x, 1]))+2f((Subscript[x, 2]))+4f((Subscript[x, 3]))+f((Subscript[x, 4])))
so
	(∫_0^1)(e^(-(x^2)))d x≈(0.25/3)((e^(-(((0))^2)))+4(e^(-(((0.25))^2)))+2(e^(-(((0.5))^2)))+4(e^(-(((0.75))^2)))+(e^(-(((1))^2))))
	≈0.746855

We can use the error estimate to determine the error in our approximation. Using K=12,a=0,b=1,n=4, Simpson’s error estimate is
	|exact-0.746855|≤(((K(b-a))^5)/(180(n^4)))=((12(((1-0))^5))/(180×(4^4)))≈0.00026
We know our approximation to the integral is correct to the first 3 decimal places.

Question

The definite integral (∫_0^1)sin((Sqrt[1+(x^4)]))d x is an integral for which we do not know an antiderivative. Approximate its value with error less than or equal to 0.0001 using Simpson’s rule. (Remember to use an even number n=2m of subintervals.)

Note: If your call to simpsonsRule invokes n subintervals for midpoint and trapezoidal rules, then you would need to call simpsonsRule with m=n/2, rather than n. Remember: you are using your superpower to create the rule: you need to live with your work!

simpsonsRule[Sin[Sqrt[1+x^4]],{x,0,1,4}] (* If I found 8 for n, my implementation of simpsonsRule would demand I use 4 for the call.... *)
NIntegrate[Sin[Sqrt[1+x^4]],{x,0,1}]