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

Riemann sums

Question

Approximate the area between the x-axis and y=(e^(-(x^2))) for 0≤x≤2 using four rectangles:

Plot[E^(-x^2),{x,0,2},PlotRange->{0,1}]
NIntegrate[E^(-x^2),{x,0,2}]

There are several different (common) rectangle rules:

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
Of these three, the Midpoint method is the most accurate.

Trapezoid rule

The trapezoid 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. Let Δ x=((b-a)/n) and (Subscript[x, k])=a+k Δ x.
	(∫_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])))
	The Trapezoidal method is about half as accurate as the Midpoint method: so Midpoint is still the winner (in general).

Example

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

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 approximation? Find the absolute error |exact-approximation|.

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

Questions

Is the trapezoid rule the same as the midpoint rule?

Midpoint rule:
(∫_a^b)f(x) d 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

Trapezoid rule: 
(∫_a^b)f(x) d x≈(1/2)Δ x(f((Subscript[x, 0]))+2f((Subscript[x, 1]))+…+2f((Subscript[x, n-1]))+f((Subscript[x, n])))

Use the trapezoid rule to estimate (∫_0^2)(Sqrt[1+(x^3)])d x with 5 trapezoids.

Δ x=((2-0)/5)=0.4
(Subscript[x, 0])=0
(Subscript[x, 1])=0.4
(Subscript[x, 2])=0.8
(Subscript[x, 3])=1.2
(Subscript[x, 4])=1.6
(Subscript[x, 5])=2.0

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|.

Questions

Use the trapezoid rule with n=5 to estimate (∫_0^1)(x^2)d x.

Compute the exact value and determine the absolute error in using this approximation.

Δ x=((1-0)/5)=0.2; (Subscript[x, 0])=0,(Subscript[x, 1])=0.2,(Subscript[x, 2])=0.4,(Subscript[x, 3])=0.6,(Subscript[x, 4])=0.8,(Subscript[x, 5])=1

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

trapezoid[(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!

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.

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.