Numerical Integration

MAT 229, Spring 2021

Week 6

Supporting materials

If you wish to get a different perspective on the notes below, try either of the following textbook sections.

Stewart’s Calculus
	Section 7.7: Approximate integration

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

Rectangles

Consider the integral (∫_a^b)f(x)d x. To approximate the area under the curve y=f(x) with rectangles of equal width, two things are needed:

riemannSumPlot[exp_,{x_,xmin_,xmax_,n_},opts___]:=Module[{func,dx,type,offset,val,var,plot,rectangles,intVal},
type = RiemannType/.{opts}/.RiemannType->"MidPoint";
If[type=="RightPoint",offset = 0.0,
If[type=="MidPoint",offset = 0.5,
If[type=="LeftPoint",offset = 1.0,
offset=type
]
]
];
dx = N[(xmax-xmin)/n];
func[t_]:=exp/.x->t;
var=Table[(xmin+dx(offset+k)),{k,0,n-1}];
val=Map[func,var];
intVal=Total[val]*dx;
plot=Plot[exp,{x,xmin,xmax}];
rectangles=Graphics[Join[{Opacity[0.2],Gray},Table[Rectangle[{xmin+k dx,0},{xmin+(k+1)dx,val[[k+1]]}],{k,0,n-1}]]];
Show[rectangles,plot,FilterRules[{opts},Options[Graphics]],Axes->True,PlotLabel->intVal]
]

Show[Plot[2+x Cos[x],{x,-2,3}],riemannSumPlot[2+x Cos[x],{x,-1,2,6},RiemannType->0.7,AspectRatio->1],PlotRange->{0,3},AxesOrigin->{0,0},PlotLabel->Integrate[f[x],{x,-1,2}]]

The number of rectangles.

Where each rectangle should touch the curve, i.e. which value of x to use in the function evaluation (rectangle height) in each subinterval.

Rectangle base points

The entire width of this region is b-a. If we decide on n rectangles, the width of each rectangle is
	Δ  x=((b-a)/n)
The locations of the base of each rectangle along the x-axis:

The first rectangle on the left end starts at the lower limit of the integral and goes for one width: (Subscript[x, 0])=a,(Subscript[x, 1])=a+Δ x

The second rectangle begins where the first one ended and goes for another width: (Subscript[x, 1])=a+Δ x,(Subscript[x, 2])=a+2Δ x

The third rectangle begins where the second one ended and goes for another width: (Subscript[x, 2])=a+2Δ x,(Subscript[x, 3])=a+3Δ x

Continue the process until you reach the (n^th) rectangle: (Subscript[x, n-1])=a+(n-1)Δ x,(Subscript[x, n])=a+n Δ x=b (Why is the last one b?)

Rectangles: heights given by function values, chosen from the sub-intervals

Three standard ideas about where to sample the function lead to three different rectangle rules:  left endpoint, the midpoint, and the right endpoint.

GraphicsRow[{Show[Plot[2+x Cos[x],{x,-2,3}],riemannSumPlot[2+x Cos[x],{x,-1,2,6},RiemannType->0.0,AspectRatio->1],PlotRange->{-1,3},AxesOrigin->{0,0},PlotLabel->"Left Endpoint",Ticks->{Range[-2,3],Range[-1,3]}],Show[Plot[2+x Cos[x],{x,-2,3}],riemannSumPlot[2+x Cos[x],{x,-1,2,6},RiemannType->0.5,AspectRatio->1],PlotRange->{-1,3},AxesOrigin->{0,0},PlotLabel->"Midpoint",Ticks->{Range[-2,3],Range[-1,3]}],Show[Plot[2+x Cos[x],{x,-2,3}],riemannSumPlot[2+x Cos[x],{x,-1,2,6},RiemannType->1,AspectRatio->1],PlotRange->{-1,3},AxesOrigin->{0,0},PlotLabel->"Right Endpoint",Ticks->{Range[-2,3],Range[-1,3]}]}]

So the approximate area of a method will be a sum of little rectangular areas that look like f( (SuperStar[x]))Δ x, height*width, one for each of the n subintervals.  (SuperStar[x])  is a particular choice of  x within each subinterval (e.g. left endpoint, midpoint, right endpoint -- but other choices are possible, although we won’t get into them).

Question

Do you see that the above three rules translate to the following formulas? One thing we typically do is “distribute out” a Δ x, since the width is the same for each rectangle.

Left endpoint rule: (∫_a^b)f(x) d x≈Δ x((∑_(k=1)^n)f(a+(k-1)Δ x))=Δ x(f(a)+f(a+Δ x)+…+f(a+(n-1)Δ x))

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

Right endpoint rule: (∫_a^b)f(x) d x≈Δ x((∑_(k=1)^n)f(a+k Δ x))=Δ x(f(a+Δ x)+…+f(a+(n-1)Δ x)+f(a+n Δ x))

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. We can either 
a. calculate Left and Right separately, and then average those two; or we can 
b. average their formulas: 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])))
	=(1/2)Δ x(f(a)+2(Underoverscript[Σ, (k=1), (n-1)])f(a+k Δ x)+f(b))

A trapezoid has the following shape and its area is ((((Subscript[h, 1])+(Subscript[h, 2]))w)/2)

Graphics[{Arrowheads[0.1],Gray,Polygon[{{0,0},{1,0},{1,3},{0,2}}],Text[(Subscript[h, 1]),{-0.2,1}],Arrow[{{-0.2,1.1},{-0.2,2}}],Arrow[{{-0.2,0.8},{-0.2,0}}],Text[(Subscript[h, 2]),{1.2,1.5}],Arrow[{{1.2,1.7},{1.2,3}}],Arrow[{{1.2,1.3},{1.2,0}}],Text[w,{0.5,-.2}],Arrow[{{0.3,-0.2},{0.0,-0.2}}],Arrow[{{0.7,-0.2},{1,-0.2}}]},BaseStyle->FontSize->14]

The term “trapezoid rule” is used for this rule because the above average of the left and right endpoint rules is equivalent to summing the areas of trapezoids -- trapezoids which I hope you recognize as the average of two rectangles!

trapezoidRulePlot[exp_,{x_,xmin_,xmax_,n_},opts___]:=Module[{func,dx,val,var,plot,trapezoids,intVal},
dx = N[(xmax-xmin)/n];
func[t_]:=exp/.x->t;
var=Table[(xmin+dx k),{k,0,n}];
val=Map[func,var];
intVal=dx/3(val[[1]]+val[[-1]]+Sum[val[[k]],{k,2,n}]);
plot=Plot[exp,{x,xmin,xmax}];
trapezoids=Graphics[Join[{Opacity[0.2],Lighter[Gray,.2],EdgeForm[Gray]},Table[Polygon[{{var[[k]],val[[k]]},{var[[k]],0},{var[[k+1]],0},{var[[k+1]],val[[k+1]]}}],{k,1,n}]]];
Show[trapezoids,plot,FilterRules[{opts},Options[Graphics]],Axes->True,PlotLabel->intVal]
]

Show[Plot[2+x Cos[x],{x,-2,3}],trapezoidRulePlot[2+x Cos[x],{x,-1,2,6},AspectRatio->1],PlotRange->{0,3},AxesOrigin->{0,0},PlotLabel->Integrate[f[x],{x,-1,2}]]

Example: you don’t need a function if you have regularly spaced data.

A functional expression isn’t necessary if you have regularly spaced data:

	(({{t(hours), 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13}, {p(t)=proportion in rigor mortis , (2/123), (16/123), (47/123), (61/123), (81/123), (92/123), (99/123), (103/123), (110/123), (111/123), (111/123), (113/123)}}))

So while we may think “p(7)=(92/123)”, there’s really nothing to compute with. We just use the data. To estimate the average value of the underlying proportion p(t), 2≤t≤13, estimate the interval (1/(13-2))(∫_2^13)p(t)d t with the trapezoid rule on the data using n=11.
And while it might be better to speak of t rather than x, and “Δ t” rather than Δ x, we’ll continue to use x and Δ x:

Δ x=((13-2)/11)=1

(Subscript[x, 0])=2,(Subscript[x, 1])=3,(Subscript[x, 2])=4,…,(Subscript[x, 11])=13

(1/11)(∫_2^13)p(t)d t≈ (1/11)((1/2)Δ x(p((Subscript[x, 0]))+2p((Subscript[x, 1]))+…+2p((Subscript[x, 10]))+p((Subscript[x, 11]))))
		=(1/2)(1)(({{(2/123)+2((16/123), +(47/123), +(61/123)+, (81/123)+, (92/123)+, (99/123)+, (103/123)+, (110/123)+, (111/123)+, (111/123)), + (113/123)}}))
				=(1/22)((1/123))(2+32+94+122+162+184+198+206+220+222+222+113)
				=(1/246)(1777)≈0.6566888396156689`

Technology

For values of n other than small ones, using these numerical techniques by hand is tedious. However, for computing devices it is a piece of cake. Two things make the process easier to implement.

A simple way to represent a function, as in the function notation f(x)=x (e^(-2x)).

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

A simple way to implement the summation notation: (Underoverscript[Σ, (k=1), n])g(k)=g(1)+g(2)+g(3)+…+g(n-1)+g(n).

a=0.0
b=1.0
n=10
dx=(b-a)/n (* The width of a subinterval. *)
g[k_]:=f[a+(k-1)dx]dx (* Areas of kth rectangle over kth subintervals, Left Rectangle Rule *)
Sum[g[k],{k,1,10}]  (* Add up the areas, and get a good LRR approximation to the integral: *)
Integrate[f[x],{x,a,b}]

Example

Consider the integral (∫_2^6)(Sqrt[(x^3)-5])d x.

Use the midpoint rule with n=50 to approximate it. (∫_a^b)f(x) d x≈Δ x((∑_(k=1)^n)f(a+(k-(1/2))Δ x))

f[x_]:=Sqrt[x^3-5]

a=2.0
b=6.0
n=50
dx=(b-a)/n

dx*Sum[f[a+(k-1/2)*dx],{k,1,n}]

Estimate the error using the midpoint rule error estimate, absolute error≤(((K(b-a))^3)/(24(n^2))) where K≥|f''(x)| for a≤x≤b.

Plot[Abs[f''[x]],{x,a,b},PlotRange->{0,4}]

From this we can take K=4. (Could even do 3.5, but just need to get a value at least as big. It’s always better to err on the safe side.)

K2=4;
K2(((b-a))^3)/(24*n^2)

We know our estimate of 31.4195 is accurate to a couple of decimal places.

Question

Using the trapezoid rule, estimate (∫_1^4)(x^(x/2))d x with error less than 0.001.

The error estimate for the trapezoid rule is absolute error≤(((K(b-a))^3)/(12(n^2))) where K≥|f''(x)| for a≤x≤b. Define the function and find a value for K.

g[x_]:=x^(x/2)

a=1.0;
b=4.0;
Plot[Abs[g''[x]],{x,a,b}]

We can use the right end point height for K.

K2=Abs[g''[b]]

Just to be on the safe side, let’s round up to 25. Now solve ((25(((4-1))^3))/(12(n^2)))=0.001 for n.

Clear[n]
K2=25;
NSolve[K2(b-a)^3/(12n^2)==0.001,n]

Use the positive value and round up to a whole number: n=238. The final estimate is 
	(∫_a^b)f(x) d x≈(1/2)Δ x(f(a)+2(Underoverscript[Σ, (k=1), (n-1)])f(a+k Δ x)+f(b))

n=238;
dx=(b-a)/n

0.5dx(g[a]+2Sum[g[a+k*dx],{k,1,n-1}]+g[b])

Our estimate is (∫_1^4)(x^(x/2))d x≈14.1914 and we know it is accurate to the number of decimal places displayed.

(Work the problem, then compare your results with my results by clicking on the ⟩⟩ button above.)

Simpson’s rule

The error estimates for the trapezoid rule and the midpoint rule are similar.

Trapezoid rule error estimate: absolute error≤(((K(b-a))^3)/(12(n^2)))

Midpoint rule error estimate: absolute error≤(((K(b-a))^3)/(24(n^2)))

The midpoint rule looks to be about twice as accurate as the trapezoid rule using the same n. Simpson's rule for estimating (∫_a^b)f(x) d x  is the weighted average of the midpoint rule (counted twice) and the trapezoid rule (counted once), but computed with m=2n subintervals.

(Subscript[Simpsons, m])=((2(Subscript[Midpoint, n])+(Subscript[Trapezoid, n]))/3)

This is the easy way to compute Simpson’s rule! However, you could perform this weighted average on the rules themselves:

Let Δ x=((b-a)/m) and (Subscript[x, k])=a+k Δ x.
	(∫_a^b)f(x) d x≈(1/3)Δ x(f((Subscript[x, 0]))+4f((Subscript[x, 1]))+2f((Subscript[x, 2]))+4f((Subscript[x, 3]))+2f((Subscript[x, 4]))+…+2f((Subscript[x, 2n-2]))+4f((Subscript[x, 2n-1]))+f((Subscript[x, 2n])))
	
Note that except for the first and last terms in this sum

the odd indexed terms have a coefficient of 4: index=2k-1

the even indexed terms have a coefficient of 2: index=2k

In summation notation, Simpson’s rule is (∫_a^b)f(x) d x≈(1/3)Δ x (f(a)+4 (Underoverscript[Σ, (k=1), n])f((Subscript[x, 2k-1]))+2(Underoverscript[Σ, (k=1), (n-1)])f((Subscript[x, 2k])) + f(b))

Example

Consider the integral (∫_2^6)(Sqrt[(x^3)-5])d x.

Use Simpson’s rule with m=2n=100 (so that n=50) to approximate it. (∫_a^b)f(x) d x≈(1/3)Δ x (f(a)+4 (Underoverscript[Σ, (k=1), n])f(a+(2k-1)Δ x)+2(Underoverscript[Σ, (k=1), (n-1)])f(a+2k Δ x)+f(b))

f[x_]:=Sqrt[x^3-5]

a=2.0
b=6.0
n=50
m=2*n
dx=(b-a)/m

(1/3)dx*(f[a]+4Sum[f[a+(2k-1)dx],{k,1,n}]+2Sum[f[a+2k*dx],{k,1,n-1}]+f[b])

Estimate the error using (Subscript[S, m]) ,absolute error≤(((K(b-a))^5)/(180(m^4))) where K≥|(f^((4)))(x)| for a≤x≤b.

Plot[Abs[f''''[x]],{x,a,b},PlotRange->{0,210}]

It looks like the highest point is at x=a.

K4=Abs[f''''[a]]

From this we can take K=201.

K4=201
K4(((b-a))^5)/(180*m^4)

We know our estimate of 31.4196 is accurate to at least four decimal places.

Question

Using Simpson’s rule, estimate (∫_1^4)(x^(x/2))d x with error less than 0.001.

The error estimate for Simpson’s rule (Subscript[S, m]) absolute error≤(((K(b-a))^5)/(180(m^4)))  where K≥|(f^((4)))(x)| for a≤x≤b. Define the function and find a value for K. Notice that we are thinking of m=2n, an even number. Find the corresponding value of m to assure an error less than 0.001.

g[x_]:=x^(x/2)

a=1.0;
b=4.0;
Plot[Abs[g''''[x]],{x,a,b}]

We can use the right end point height for K.

K4=Abs[g''''[b]]

Just to be on the safe side, let’s round up to 49. Now solve ((49(((4-1))^5))/(180(m^4)))=0.001 for m.

Clear[m]
K4=49
NSolve[K4(b-a)^5/(180m^4)==0.001,m,Reals]

Use the positive value and round up to the next even whole number: m=18=2n, n=9. The final estimate is 
	(∫_a^b)f(x) d x≈(1/3)Δ x (f(a)+4 (Underoverscript[Σ, (k=1), n])f(a+(2k-1)Δ x)+2(Underoverscript[Σ, (k=1), (n-1)])f(a+2k Δ x)+f(b))

So we could compute (Subscript[T, 9]) and (Subscript[M, 9]), and average them; or we could proceed directly:

n=9
m=2n
dx=(b-a)/m

1/3dx(g[a]+4Sum[g[a+(2k-1)*dx],{k,1,n}]+2Sum[g[a+2k*dx],{k,1,n-1}]+g[b])

Our estimate is (∫_1^4)(x^(x/2))d x≈14.191 and we know it is accurate to the number of decimal places displayed.

(Work the problem, then compare your results with my results by clicking on the ⟩⟩ button above.)

Perfect integration from Numerical Schemes

Because of the form of the error bounds for the Trapezoidal, Midpoint , and Simpson’s rules, you should realize that the Trapezoidal and Midpoint rules get integrals of linear functions exactly right; and that Simpson’s rule gets integrals of cubic polynomials exactly right. No matter how many subintervals you take!

Let’s see.... The fewest subintervals Simpson’s can use is 2. So take m=2 (then n=1 for the Trapezoidal and Midpoint rules).

Example

Evaluate the definite integral (∫_0^3)(16+(x^2)-(x^3))d x

a=0.0
b=3.0
n=1
dx=(b-a)/n
f[x_]:=16+x^2-x^3
Plot[f[x],{x,0,3}]

Integrate[f[x],{x,0,3}]
true=N[%]

lrrN=dx*Sum[f[a+(k-1)dx],{k,1,n}]
lrrN-true (* dramatic overestimate)

rrrN=Sum[dx*f[a+(k-0) dx],{k,1,n}]
rrrN-true  (* dramatic underestimate)

midN=dx*Sum[f[a+(k-1/2) dx],{k,1,n}]
midN-true

trapN=1/2(lrrN+rrrN)
trapN-true

Error estimate

For the trapezoid rule the absolute error is less than or equal to
	(((K(b-a))^3)/(12(n^2)))

Plot[Abs[f''[x]],{x,a,b}]

K2=Abs[f''[b]]

trapErr=K2(b-a)^3/(12*n^2)

midErr=K2(b-a)^3/(24*n^2)

Estimate and Error estimate for Simpson’s Rule:

S2N=(2*midN+trapN)/3
true
S2N-true

Perfect! From a bunch of crummy estimates, Simpson’s got it exactly right! How?

For Simpson’s rule the absolute error is less than or equal to
	(((K(b-a))^5)/(180(n^4)))

Plot[Abs[f''''[x]],{x,0,3}] (* identically zero! *)

K4=0 
(* Since the fourth derivative is identically zero! *)

K4(b-a)^5/(180*(2n)^4) 
(* = 0; Simpson's gets it exactly right! *)