Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
187 views

#Numerical Integration

Some antiderivatives are difficult (or impossible) to compute. In such cases, numerical approximation of a definite integral may be a better (or the only) option.

In Calculus 1, we learned about the "numerical_integral" command in Sage. We are not going to use this command in this lab. Instead, we will explore various approximation techniques.

##Riemann Sums

The simplest numerical approximation comes from the definition of the definite integral as a limit of Riemann sums. Each Riemann sum is an approximation of the definite integral using rectangles (with one side on the x-axis), and we can make the approximation better by increasing the number of rectangles.

Our goal is to approximate abf(x)dx\displaystyle\int_a^bf(x)\,dx.

We will use nn rectangles of equal width, which means the width of each rectangle is Δx=ban\Delta x=\frac{b-a}{n}.

The height of the rectangle will be given by the value of the function ff at some point in the base of that rectangle.

If we choose the left endpoint of each rectangle, we will find the left Riemann sum. If we choose the right endpoint of each rectangle, we will find the right Riemann sum.

The endpoints of the rectangles are a, a+Δx, a+2Δx, a+3Δx,..., a+nΔx=ba,\ a+\Delta x,\ a+2\Delta x,\ a+3\Delta x,... ,\ a+n\Delta x=b.

So the left Riemann sum is LS=f(a)Δx+f(a+Δx)Δx++f(a+(n1)Δx)Δx=i=0n1(f(a+iΔx)Δx)\displaystyle LS=f(a)\cdot\Delta x+f(a+\Delta x)\cdot\Delta x+\cdots+f(a+(n-1)\Delta x)\cdot\Delta x=\sum_{i=0}^{n-1}\left(f(a+i\Delta x)\cdot\Delta x\right).

Similarly, the right Riemann sum is RS=f(a+Δx)Δx+f(a+2Δx)Δx++f(b)Δx=i=1n(f(a+iΔx)Δx)\displaystyle RS=f(a+\Delta x)\cdot\Delta x+f(a+2\Delta x)\cdot\Delta x+\cdots+f(b)\cdot\Delta x=\sum_{i=1}^{n}\left(f(a+i\Delta x)\cdot\Delta x\right).

The only difference in these formulas is the index of summation.

###Example 1

Consider 012x2xdx\displaystyle\int_0^12x^2-x\, dx.

Here are pictures with n=5n=5 rectangles.

f(x)=2*x^2-x #function to be integrated a=0 #lower limit of integration b=1 #upper limit of integration n=5 #number of subintervals dx=(b-a)/n #delta x p=plot(f,(a,b),color='black');q=p for i in [0..n-1]: #adds rectangles for left sum p=p+polygon([(a+i*dx,0),(a+i*dx,f(a+i*dx)),(a+(i+1)*dx,f(a+i*dx)),(a+(i+1)*dx,0)],alpha=.5)+polygon([(a+i*dx,0),(a+i*dx,f(a+i*dx)),(a+(i+1)*dx,f(a+i*dx)),(a+(i+1)*dx,0)],fill=False,color='gray') p.show(title='Left Riemann Sum') for i in [0..n-1]: #adds rectangles for right sum q=q+polygon([(a+i*dx,0),(a+i*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,0)],alpha=.5)+polygon([(a+i*dx,0),(a+i*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,0)],fill=False,color='gray') q.show(title='Right Riemann Sum')

The formulas below compute the left and right Riemann sums.

f(x)=2*x^2-x #function to be integrated a=0 #lower limit of integration b=1 #upper limit of integration n=5 #number of subintervals dx=(b-a)/n #delta x i=var('i') LS=sum(f(a+i*dx)*dx,i,0,n-1) #calculates left sum print 'The left Riemann sum is',LS,'=',LS.n(digits=5) RS=sum(f(a+i*dx)*dx,i,1,n) #calculates right sum print 'The right Riemann sum is',RS,'=',RS.n(digits=5)
The left Riemann sum is 2/25 = 0.080000 The right Riemann sum is 7/25 = 0.28000

For the sake of comparison, the exact value is 160.1666667\frac{1}{6}\approx0.1666667.

If you want to graph or compute any other left or right sums, you may copy and paste the above formulas and change ff, aa, bb, and nn. Leave everything else!

Let's change the number of rectangles to 10 and see how our approximation improves.

f(x)=2*x^2-x a=0 b=1 n=10 #Just changed this to 10 dx=(b-a)/n p=plot(f,(a,b),color='black');q=p for i in [0..n-1]: p=p+polygon([(a+i*dx,0),(a+i*dx,f(a+i*dx)),(a+(i+1)*dx,f(a+i*dx)),(a+(i+1)*dx,0)],alpha=.5)+polygon([(a+i*dx,0),(a+i*dx,f(a+i*dx)),(a+(i+1)*dx,f(a+i*dx)),(a+(i+1)*dx,0)],fill=False,color='gray') p.show(title='Left Riemann Sum') for i in [0..n-1]: q=q+polygon([(a+i*dx,0),(a+i*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,0)],alpha=.5)+polygon([(a+i*dx,0),(a+i*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,0)],fill=False,color='gray') q.show(title='Right Riemann Sum') print '' i=var('i') LS=sum(f(a+i*dx)*dx,i,0,n-1) #calculates left sum print 'The left Riemann sum is',LS,'=',LS.n(digits=5) RS=sum(f(a+i*dx)*dx,i,1,n) #calculates right sum print 'The right Riemann sum is',RS,'=',RS.n(digits=5)
The left Riemann sum is 3/25 = 0.12000 The right Riemann sum is 11/50 = 0.22000

Remember that the exact answer is 16\frac{1}{6}. Both of these are closer to the correct answer than before, but they are still not very close.

Of course, we can increase our accuracy by using larger values of nn, but it would be very nice if there were more accurate approximation methods. Fortunately, there are.

##Midpoint and Trapezoidal Rules

The Midpoint Rule is another Riemann sum approximation, but instead of using the left or right endpoints, we will use the midpoint of each rectangle.

Using the midpoint of each rectangle as a height often leads to better estimates with a similar number of rectanges compared to the left or right points. Notice that using right endpoints an increasing function will be consistantly overestimated and left end points underestimated. With midpoints these over and underestimates tend to cancel giving a better overall approximation.

###Example 2

Consider 012x2xdx\displaystyle\int_0^1 2x^2-x\, dx.

Here is a plot and calculation for the Midpoint Rule using n=5n=5.

f(x)=2*x^2-x #function to be integrated a=0 #lower limit of integration b=1 #upper limit of integration n=5 #number of subintervals dx=(b-a)/n #delta x m=plot(f,(a,b),color='black') for i in [0..n-1]: #sets up the plot m=m+polygon([(a+i*dx,0),(a+i*dx,f(a+i*dx+dx/2)),(a+(i+1)*dx,f(a+i*dx+dx/2)),(a+(i+1)*dx,0)],alpha=.5)+polygon([(a+i*dx,0),(a+i*dx,f(a+i*dx+dx/2)),(a+(i+1)*dx,f(a+i*dx+dx/2)),(a+(i+1)*dx,0)],fill=False,color='gray') m.show(title='Midpoint Rule') print '' i=var('i') M=sum(f(a+i*dx+dx/2)*dx,i,0,n-1) #performs the calculation print 'The Midpoint Rule gives',M,'=',M.n(digits=5)
The Midpoint Rule gives 4/25 = 0.16000
%md The exact value is $\frac{1}{6}$, and the Midpoint Rule gives us the correct answer to two decimal places with only 5 rectangles. Recall that for $n=5$, the left sum was 0.12 and the right sum was 0.22. The Midpoint Rule has done a good job of balancing out the under- and over-estimates. However, if our function was not so nicely behaved (here the function is similar on bothsizes of the midpont) the midpoint rule may not benefit from the cancellation effect.

The exact value is 16\frac{1}{6}, and the Midpoint Rule gives us the correct answer to two decimal places with only 5 rectangles.

Recall that for n=5n=5, the left sum was 0.12 and the right sum was 0.22. The Midpoint Rule has done a good job of balancing out the under- and over-estimates. However, if our function was not so nicely behaved (here the function is similar on bothsizes of the midpont) the midpoint rule may not benefit from the cancellation effect.

So far we have been approximating our function using a constant (horizontal line = degree 0 polynomial) on each subinterval.

We may be able to improve our approximation if we use a degree 1 polynomial instead (i.e., a non-horizontal line). For each subinterval, we’ll use the secant line based on the left and right endpoints. Of course, the resulting shape is not a rectangle but a trapezoid. Thus, this approach is called the Trapezoidal Rule.

Recall that the area of a trapezoid is b1+b22h\displaystyle \frac{b_{1}+b_{2}}{2}h where b1b_{1} and b2b_{2} are the lengths of the bases and hh the height. In this case, the trapezoid is sitting on its side, so the bases are actually vertical and the height is Δx\Delta x.

###Example 3

Consider 012x2xdx\displaystyle \int_{0}^{1}2x^{2}-xdx

Here is the plot and calculation for the Trapezoidal Rule using n=5n=5.

f(x)=2*x^2-x #function to be integrated a=0 #lower limit of integration b=1 #upper limit of integration n=5 #number of subintervals dx=(b-a)/n #delta x t=plot(f,(a,b),color='black') for i in [0..n-1]: #sets up the plot t=t+polygon([(a+i*dx,0),(a+i*dx,f(a+i*dx)),(a+(i+1)*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,0)],alpha=.5)+polygon([(a+i*dx,0),(a+i*dx,f(a+i*dx)),(a+(i+1)*dx,f(a+(i+1)*dx)),(a+(i+1)*dx,0)],fill=False,color='gray') t.show(title='Trapezoidal Rule') print '' i=var('i') T=sum((f(a+i*dx)+f(a+(i+1)*dx))*dx/2,i,0,n-1) #performs the calculation print 'The Trapezoidal Rule gives',T,'=',T.n(digits=5)
The Trapezoidal Rule gives 9/50 = 0.18000

Notice that this approximation is slightly worse than what we got from the Midpoint Rule, but it is much better than either the left or right sum.

In fact, numerically the Trapezoidal Rule simply gives the average of the left and right sums!

In general (but not always), the Midpoint Rule will give a better approximation. If the function ff is either concave up or concave down on the entire subinterval, then the linear approximation from the Trapezoidal Rule (the secant line) will be entirely above or below the curve. On the other hand, the horizontal line segment from the Midpoint Rule will be above the curve on part of the interval and below the curve on part of the interval, so the errors will cancel out.

Despite this numerical disappointment, the idea of increasing the degree of the approximating polynomial was sound. The next step is to use a degree 2 polynomial (quadratic, parabola) to approximate ff on each subinterval. This is called Simpson's Rule.

##Simpson's Rule

Instead of using a line to approximate our curve on each subinterval, we will use a parabola. Since a parabola is usually closer to the curve than a line, this should give us a better approximation (of course, this is not true if our original curve is a line, but in that case, why are we approximating the integral?).

A line is determined by 2 points, but it takes 3 points to determine a parabola.

The three points that are normally used are the left endpoints of three consecutive subintervals. This choice forces us to use an even number of subintervals (i.e., nn must be even).

I won't take you through all the calculations (see the textbook, pages 463-465).

###Example 4

Consider 012x2xdx\displaystyle\int_0^1 2x^2-x\, dx.

Here is the result for Simpson's Rule using n=4n=4. (I'll skip the graph. If you really want to see one, there is an example at the end of this worksheet.)

f(x)=2*x^2-x #function to be integrated a=0 #lower limit of integration b=1 #upper limit of integration n=4 #number of subintervals (must be even) dx=(b-a)/n #delta x n2=int(n/2) coeffs = [4,2]*n2 coeffs = [1] +coeffs[:n-1]+[1] S=(dx/3)*sum([coeffs[i]*f(a+dx*i) for i in [0..n]]) print 'Simpson\'s Rule gives',S,'=',S.n(digits=5)
Simpson's Rule gives 1/6 = 0.16667

Notice we obtain an exact answer utilizing Simpson's Rule. This should not be too surprizing as our function was a parabola in the first place.

︠ec3d9033-b2d2-41d7-a9ce-278fdab62d93i︠ %md It is interesting to note that numerically the answer from Simpson's Rule is a weighted average of the answers from the Trapezoidal Rule and the Midpoint Rule (with half the number of rectangles). In particular, if $S_{2n}$ is the approximation from Simpson's Rule with $2n$ subintervals, $T_n$ is the approximation from the Trapezoidal Rule with $n$ trapezoids, and $M_n$ is the approximation from the Midpoint Rule with $n$ rectangles, then $$S_{2n}=\frac{T_n+2M_n}{3}$$

It is interesting to note that numerically the answer from Simpson's Rule is a weighted average of the answers from the Trapezoidal Rule and the Midpoint Rule (with half the number of rectangles). In particular, if S2nS_{2n} is the approximation from Simpson's Rule with 2n2n subintervals, TnT_n is the approximation from the Trapezoidal Rule with nn trapezoids, and MnM_n is the approximation from the Midpoint Rule with nn rectangles, then S2n=Tn+2Mn3S_{2n}=\frac{T_n+2M_n}{3}

###Example 5 (and a one-stop shop for cut-and-paste)

Here is one example that puts all the rules together. (I'll skip the graphs again.).

Let's approximate 14x34x2+6x9dx\displaystyle\int_1^4 x^3-4x^2+6x-9\, dx using n=10n=10, n=50n=50, and n=100n=100.

For comparison, the exact answer is 94=2.25-\frac{9}{4}=-2.25.

f(x)=x^3-4*x^2+6*x-9 #function to be integrated a=1 #lower limit of integration b=4 #upper limit of integration n=10 #number of subintervals displaydigits=5 #number of digits to display in final answer dx=(b-a)/n i=var('i') LS=sum(f(a+i*dx)*dx,i,0,n-1) #calculates left sum print 'The left Riemann sum is',LS.n(digits=displaydigits) RS=sum(f(a+i*dx)*dx,i,1,n) #calculates right sum print 'The right Riemann sum is',RS.n(digits=displaydigits) M=sum(f(a+i*dx+dx/2)*dx,i,0,n-1) print 'The Midpoint Rule gives',M.n(digits=displaydigits) T=sum((f(a+i*dx)+f(a+(i+1)*dx))*dx/2,i,0,n-1) print 'The Trapezoidal Rule gives',T.n(digits=displaydigits) n2=int(n/2) coeffs = [4,2]*n2 coeffs = [1] +coeffs[:n-1]+[1] S=(dx/3)*sum([coeffs[i]*f(a+dx*i) for i in [0..n]]) print 'Simpson\'s Rule gives',S.n(digits=displaydigits)
The left Riemann sum is -5.2425 The right Riemann sum is 1.0575 The Midpoint Rule gives -2.3288 The Trapezoidal Rule gives -2.0925 Simpson's Rule gives -2.2500

To summarize: L10=5.2425, R10=1.0575, M10=2.3288, T10=2.0925, S10=2.25L_{10}=-5.2425,\ R_{10}=1.0575,\ M_{10}=-2.3288,\ T_{10}=-2.0925,\ S_{10}=-2.25, and the exact value is 2.25-2.25.

Now let's try n=50n=50:

n=50 #number of subintervals displaydigits=5 #number of digits to display in final answer dx=(b-a)/n i=var('i') LS=sum(f(a+i*dx)*dx,i,0,n-1) #calculates left sum print 'The left Riemann sum is',LS.n(digits=displaydigits) RS=sum(f(a+i*dx)*dx,i,1,n) #calculates right sum print 'The right Riemann sum is',RS.n(digits=displaydigits) M=sum(f(a+i*dx+dx/2)*dx,i,0,n-1) print 'The Midpoint Rule gives',M.n(digits=displaydigits) T=sum((f(a+i*dx)+f(a+(i+1)*dx))*dx/2,i,0,n-1) print 'The Trapezoidal Rule gives',T.n(digits=displaydigits) n2=int(n/2) coeffs = [4,2]*n2 coeffs = [1] +coeffs[:n-1]+[1] S=(dx/3)*sum([coeffs[i]*f(a+dx*i) for i in [0..n]]) print 'Simpson\'s Rule gives',S.n(digits=displaydigits)
The left Riemann sum is -2.8737 The right Riemann sum is -1.6137 The Midpoint Rule gives -2.2532 The Trapezoidal Rule gives -2.2437 Simpson's Rule gives -2.2500

To summarize: L50=2.8737, R50=1.6137, M50=2.2532, T50=2.2437, S50=2.25L_{50}=-2.8737,\ R_{50}=-1.6137,\ M_{50}=-2.2532,\ T_{50}=-2.2437,\ S_{50}=-2.25, and the exact value is 2.25-2.25.

Now let's try n=100n=100:

n=100 #number of subintervals displaydigits=5 #number of digits to display in final answer dx=(b-a)/n i=var('i') LS=sum(f(a+i*dx)*dx,i,0,n-1) #calculates left sum print 'The left Riemann sum is',LS.n(digits=displaydigits) RS=sum(f(a+i*dx)*dx,i,1,n) #calculates right sum print 'The right Riemann sum is',RS.n(digits=displaydigits) M=sum(f(a+i*dx+dx/2)*dx,i,0,n-1) print 'The Midpoint Rule gives',M.n(digits=displaydigits) T=sum((f(a+i*dx)+f(a+(i+1)*dx))*dx/2,i,0,n-1) print 'The Trapezoidal Rule gives',T.n(digits=displaydigits) n2=int(n/2) coeffs = [4,2]*n2 coeffs = [1] +coeffs[:n-1]+[1] S=(dx/3)*sum([coeffs[i]*f(a+dx*i) for i in [0..n]]) print 'Simpson\'s Rule gives',S.n(digits=displaydigits)
The left Riemann sum is -2.5634 The right Riemann sum is -1.9334 The Midpoint Rule gives -2.2508 The Trapezoidal Rule gives -2.2484 Simpson's Rule gives -2.2500

To summarize: L100=2.5634, R100=1.9334, M100=2.2508, T100=2.2484, S100=2.25L_{100}=-2.5634,\ R_{100}=-1.9334,\ M_{100}=-2.2508,\ T_{100}=-2.2484,\ S_{100}=-2.25, and the exact value is 2.25-2.25.

What do we see from this example?

  1. All the approximations improve as nn increases (except for Simpson's, which gives the exact answer for any cubic polynomial).

  2. The left and right Riemann sums are not very good approximations.

  3. The Midpoint Rule and Trapezoidal Rule are better than the left and right sums.

  4. The Midpoint Rule is better than the Trapezoidal Rule.

  5. Simpson's Rule is the best.

###Example 6 - Simpson's Rule Graph

The function to be integrated in in black. The red vertical lines mark the ends of each parabola. The approximating parabolas are in blue.

I have set n=10n=10. Try changing this and see how fast the parabolas get close to the original curve.

f(x)=e^x*cos(2*x) #function to be integrated a=0 #lower limit of integration b=2*pi #upper limit of integration n=10 #number of subintervals dx=(b-a)/n v=[a] w=[f(a)] p=plot(f,(a,b),color='black') for i in [0..n-1]: v=v+[a+(i+1)*dx] w=w+[f(a+(i+1)*dx)] for i in [0,2..n-1]: A,B,C=var('A B C') eqn1=A*v[i]^2+B*v[i]+C==w[i] eqn2=A*v[i+1]^2+B*v[i+1]+C==w[i+1] eqn3=A*v[i+2]^2+B*v[i+2]+C==w[i+2] coeff=solve([eqn1,eqn2,eqn3],A,B,C) A=coeff[0][0].rhs();B=coeff[0][1].rhs();C=coeff[0][2].rhs() p+=plot(A*x^2+B*x+C,(v[i],v[i+2]),fill='axis')+line([(v[i],0),(v[i],f(v[i]))],color='red',linestyle='--') p+=line([(v[n],0),(v[n],f(v[n]))],color='red',linestyle='--');p

#Numerical Integration Assignment

###Question 1

The Golden Gate Bridge has a main span of 4,200 feet (the distance between the two towers). The main suspension cables that support the road over this span each form a parabolic shape. The length of each cable is found by Length=0212001+0.000577x2dx\mathrm{Length}=\int_0^{21}200\sqrt{1+0.000577x^2}\, dx

####Part a

Approximate the value of this integral using left and right Riemann sums, the Midpoint Rule, the Trapezoidal Rule, and Simpson's Rule using n=10n=10, n=20n=20, and n=30n=30.

####Part b

What is your best approximation for the length of each cable (round to one decimal place)? How confident are you in this approximation?

###Question 2 Consider the integral 121x2dx\int_{-1}^{\sqrt{2}}\frac{1}{x^2}\, dx

####Part a

Use the Midpoint Rule with n=10n=10 to approximate this integral.

####Part b

Approximate this integral with the Midpoint Rule using three more values of nn.

####Part c

What happens to your approximations as nn increases? What do you conclude about the integral?

###Question 3

Consider the function f(x)=x5e2x\displaystyle f(x)=x^5e^{-2x} over the interval [2,2][-2,2].

####Part a

Approximate the area under this curve using the Midpoint Rule, Trapezoidal Rule, and Simpson's Rule using n=50n=50, n=100n=100, n=200n=200, and n=400n=400.

####Part b

Rounding to one decimal place, this area is actually 363.2-363.2. For each of the three rules, which of the four values of nn give the correct answer (when rounding to one decimal place)?

###Question 4

Consider the function f(x)=1x2+1\displaystyle f(x)=\frac{1}{x^{2}+1} over the interval [1000,1000][-1000,1000]. Using our appoximation routines, what well known number is the value of 100010001x2+1dx\displaystyle \int_{-1000}^{1000}\frac{1}{x^{2}+1}dx close to? Hint: The interval width is 20002000, pick an appropriate number of approximating intervals!

︠944da2e8-b810-4023-acdc-8753604bb6c0︠