Approximating Definite Integrals restart; The restart command is used to reset the entire Maple system and clears all variables. Next, using the with statement,, we load two packages of extra commands: the student package giving the extra calculus commands shown below, and the plots package giving extra graphing commands. with(student); with(plots); Our functions and their integrals. We define two functions as Maple expressions, both of which we intend to integrate from 1 to 3. We choose the first, NiMlImZH, to be increasing and concave up on the interval, and the second, NiMlImdH, to be decreasing and concave down. f:=sqrt(4+x^4); g:=-sqrt(4+x^4)+12; We plot the functions over the interval from 0 to 3, with NiMlImZH red and NiMlImdH green. p0:=plot(0,x=0..3): p1:=plot(f,x=1..3,color=red): p2:=plot(g,x=1..3,color=green): display(p0,p1,p2); Next we integrate the functions from 1 to 3, with NiMlImZH first. int(f,x=1..3); Looks kind of crazy. We approximate this integral to 9 decimal places using the evalf command. actualf:=evalf(%); We next integrate NiMlImdH. int(g,x=1..3); actualg:=evalf(%); Approximating using the left and right rules. The reason we need to learn how to approximate integrals is because some functions, such as NiMtJSRleHBHNiMqJCklInhHIiIjIiIi, do not have an elementary antiderivative, thus denying us the use of the first fundamental theorem of calculus. In other cases, there may be an elementary antiderivative, but it may be very difficult to obtain. Now suppose we are approximating NiMtJSRpbnRHNiQtJSJmRzYjJSJ4Ry9GKTslImFHJSJiRw== . In all of the methods discussed here, we partition the interval [NiMlImFH, NiMlImJH] into NiMlIm5H equal subintervals and let NiMvKiYlJkRlbHRhRyIiIiUieEdGJiomLCYlImJHRiYlImFHISIiRiYlIm5HRiw=. The partition points are then lableled NiMvJSJhRyYlInhHNiMiIiI=, NiMmJSJ4RzYjIiIj, ..., NiMvJiUieEc2IyUibkclImJH. The left rule uses the left hand endpoint of each interval to determine the height of each of the NiMlIm5H rectangles. With NiMvJSJuRyIjNQ==, lets see what this looks like graphically by using leftbox. leftbox(f,x=1..3,10); Our approximation is the area of the 10 rectangles. The general formula is NiMtJSRzdW1HNiQqKC0lImZHNiMmJSJ4RzYjJSJpRyIiIiUmRGVsdGFHRi5GK0YuL0YtOyIiISwmJSJuR0YuRi4hIiI=. In Maple, we implement this sum with leftsum. We also approximate this answer to 9 decimal places. leftsum(f,x=1..3,10); left10:=evalf(%); In general, we define the error by error = actual value - approximation. With an increasing function, it is clear that the left rule will give an under-approximation. error_l10:=actualf-left10; Next we consider the right rule, where we use the right hand endpoint of each interval to determine the height of each of the NiMlIm5H rectangles. We again use NiMvJSJuRyIjNQ==, and use rightbox to view the situation graphically. rightbox(f,x=1..3,10); The general formula is NiMtJSRzdW1HNiQqKC0lImZHNiMmJSJ4RzYjJSJpRyIiIiUmRGVsdGFHRi5GK0YuL0YtO0YuJSJuRw==. In Maple, we implement this sum with rightsum. We again approximate this answer to 9 decimal places and find the error, which for the right rule will always be an over-approximation for an increasing function.. rightsum(f,x=1..3,10); right10:=evalf(%); error_r10:=actualf-right10; Now let's redo the above with NiMvJSJuRyIkKyI= rectangles. First , the left sum. leftbox(f,x=1..3,100); leftsum(f,x=1..3,100); left100:=evalf(%); error_l100:=actualf-left100; We see we get much less error since the area of the rectangles is much closer to the area under the curve. Next we see the same thing with the right sum. rightbox(f,x=1..3,100); rightsum(f,x=1..3,100); right100:=evalf(%); error_r100:=actualf-right100; We see that incresing the number of rectangles gives us a better approximation. Now, let's look at NiMlImdH with the left rule and NiMvJSJuRyIjNQ==. leftbox(g,x=1..3,10); It is clear that, with a decreasing fuction, the left rule will give us an over-approximation. leftsum(g,x=1..3,10); left10:=evalf(%); errorg_l10:=actualg-left10; Similarly, the right rule will give us an under-approximation. rightbox(g,x=1..3,10); rightsum(g,x=1..3,10); right10:=evalf(%); error_r10:=actualf-right10; Again, we increase the number of rectangles to NiMvJSJuRyIkKyI=, using the left rule first. leftbox(g,x=1..3,100); leftsum(g,x=1..3,100); left100:=evalf(%); error_l100:=actualg-left100; Then the right rule. rightbox(g,x=1..3,100); rightsum(g,x=1..3,100); right100:=evalf(%); error_r100:=actualg-right100; Again, more rectangles give us better approximations. Approximating using the midpoint rule. An alternate method uses the midpoint of each interval to determine the height of the rectangles to be summed. This is the midpoint rule. We view this rule graphically using NiMlImZH and middlebox. middlebox(f,x=1..3,10); The general formula is NiMtJSRzdW1HNiQqKC0lImZHNiMmJSJtRzYjJSJpRyIiIiUmRGVsdGFHRi4lInhHRi4vRi07Ri4lIm5H where NiMvJiUibUc2IyUiaUcqJiwmJiUieEdGJiIiIiZGKzYjLCZGJ0YsRiwhIiJGLEYsIiIjRjA= is the midpoint of [NiMmJSJ4RzYjLCYlImlHIiIiRighIiI=, NiMmJSJ4RzYjJSJpRw==]. In Maple, we implement this rule with middlesum. middlesum(f,x=1..3,10); mid10:=evalf(%); error_m10:=actualf-mid10; We see that we have an over-approximation here. That is because our function is concave up over the interval of integration. To see this, let's look at the midpoint rule applied to the single subinterval [2, 2.2]. middlebox(f,x=2..2.2,1); For the computations we need to do, we rewrite our function, currently a Maple expression, as a Maple function using the unapply command. ff:=unapply(f,x); This allows for easy evaluation, such as ff(2); This is in contrast to the to the following, which is necessary to find NiMtJSJmRzYjIiIj. subs(x=2,f); We use the D command to find the derivative. ffprime:=D(ff); Then we use the point-slope formula to find the tangent line to the curve at the point where the top of the rectangle meets the curve. tanline:=ff(2.1)+ffprime(2.1)*(x-2.1); We first compare the area under the rectangle (green) with the area under the tangent line (blue). p3:=plot(ff(2.1),x=2..2.2,color=green): p4:=plot(tanline,x=2..2.2,color=blue): display(p3,p4); Using congruent triangles, we see that the areas are the same. So then we compare the area under the curve (red) with the area under the tangent line. p5:=plot(f,x=2..2.2,color=red): display(p4,p5); We conclude that for a curve that is concave up, the midpoint approximation gives an under-approximation. Similarly, the midpoint approximation with a function that is concave down, such as NiMlImdH, is an over-approximation. middlebox(g,x=1..3,10); middlesum(g,x=1..3,10); mid10:=evalf(%); error_m10:=actualg-mid10; Approximating using the trapezoid rule. An alternate approach is to use the trapezoid rule, which for a given value of NiMlIm5H is the average of the left rule and righ rule, i.e., NiMvLSUldHJhcEc2IyUibkcqJiwmLSUlbGVmdEdGJiIiIi0lJnJpZ2h0R0YmRixGLCIiIyEiIg== . To see this graphically with NiMvJSJuRyIiIw==, we first define lists of our NiMlInhH and NiMtJSJmRzYjJSJ4Rw== (or NiMlInlH) values using the seq command. X:=[seq(1+1*i,i=0..2)]; Y:=[seq(ff(X[i]),i=1..3)]; We plot the graph with the trapezoids. p6:=pointplot([seq([X[i],Y[i]],i=1..3)],style=line,color=blue): p7:=pointplot({[1,0],[1,ff(1)]},style=line,color=red): p8:=pointplot({[2,0],[2,ff(2)]},style=line,color=red): p9:=pointplot({[3,0],[3,ff(3)]},style=line,color=red): display(p0,p1,p6,p7,p8,p9); We see that for a function that is concave up, the trpezoid rule will give an over-approximation. Smilarly, for a function that is concave down, the trapezoid rule gives an under-approximation. The general formula for the trapezoid rule is NiMqKCUmRGVsdGFHIiIiJSJ4R0YlIiIjISIi [NiMsKC0lImZHNiMmJSJ4RzYjIiIhIiIiKiYiIiNGKy1GJTYjJkYoNiNGK0YrRisqJkYtRistRiU2IyZGKDYjRi1GK0Yr + ... + NiMsJiomIiIjIiIiLSUiZkc2IyYlInhHNiMsJiUibkdGJkYmISIiRiZGJi1GKDYjJkYrNiNGLkYm]. It is implemented in Maple with the trapezoid command. We use the trapezoid rule with NiMvJSJuRyIjNQ==. trapezoid(f,x=1..3,10); trap10:=evalf(%); error10:=actualf-trap10; We see we get the expected over-approximation, and that it is considerably better than either the left rule or right rule, but not as good as the midpoint rule. Approximating using Simpson's rule. Since the trapezoidal error is usually about twice the midpoint error and opposite in sign, Simpson's rule, which uses a weighted average of these two, can be expected to give even more accurate approximations. The rule is NiMvLSUlc2ltcEc2IyUibkcqJiwmKiYiIiMiIiItJSRtaWRHRiZGLEYsLSUldHJhcEdGJkYsRiwiIiQhIiI= . The geometric effect of this weighted average is to use quadratic polynomials rather than straight lines to approximate curves. We illustrate this with the function NiMvLSUiRkc2IyUieEcqJkYnIiIiLSUkc2luR0YmRik= over the interval [0, 4]. We will take at first NiMvJSJuRyIiIw== to give us the three points 0, 2, 4 on the NiMlInhH-axis. Now three nonlinear points will determine a unique parabola, and the three points we use are the three points on the graph of F with NiMlInhH-coordinates 0, 2, 4. To find the formula for this parabola, we use the interp command, which takes as its arguments the lists of NiMlInhH and NiMlInlH values along with the variable chosen, here NiMlInhH. First, we directly define F as a Maple function. Notice the syntax. restart:with(student):with(plots): F:=x->x*sin(x); X:=[seq(2*i,i=0..2)]; Y:=[seq(F(X[i]),i=1..3)]; parabola:=interp(X,Y,x); Now that we have a formula for the parabola (blue), we graph it along with the function (red). p10:=plot(F(x),x=0..4,color=red): p11:=plot(parabola,x=0..4,color=blue): display(p10,p11); The parabola follows the curve pretty accurately here. Let's take the integral of the parabola and get its decimal approximation. int(parabola,x=0..4); evalf(%); Now let's approximate the integral of F by Simpson's rule by using the command simpson with NiMvJSJuRyIiIw==. Because a parabola is determined by three points covering two subintervals, NiMlIm5H must always be even when using Simpson's rule. simpson(F(x),x=0..4,2); simp2:=evalf(%); So we see that Simpson's approximation is just the integral of the parabola (with a bit of round-off error). Now let's find the integral of F and find the error in our approximation. int(F(x),x=0..4); actualF:=evalf(%); error2:=actualF-simp2; OK, this is not a real good approximation, but then our NiMlIm5H is as small as it can possibly be. Let's see what happens with NiMvJSJuRyIiKQ==. We first create a graph of the function with the four parabolas we will use. for j from 0 to 3 do X:=[seq(j+i*.5,i=0..2)]; Y:=[seq(F(X[i]),i=1..3)]; parabola:=interp(X,Y,x); pl[j]:=plot(parabola,x=j..(j+1),color=blue): od: display(p10,seq(pl[j],j=0..3)); We see that the four parabolas fit the graph of the function very well, giving hopes for a very accurate approximation. The general formula for the trapezoid rule is NiMqKCUmRGVsdGFHIiIiJSJ4R0YlIiIkISIi [NiMsKC0lImZHNiMmJSJ4RzYjIiIhIiIiKiYiIiVGKy1GJTYjJkYoNiNGK0YrRisqJiIiI0YrLUYlNiMmRig2I0YzRitGKw== +4 .NiMtJSJmRzYjJiUieEc2IyIiJA== + .. + NiMsKComIiIjIiIiLSUiZkc2IyYlInhHNiMsJiUibkdGJkYlISIiRiZGJiomIiIlRiYtRig2IyZGKzYjLCZGLkYmRiZGL0YmRiYtRig2IyZGKzYjRi5GJg==]. We now find Simpson's approximation for NiMvJSJuRyIiKQ== and the error in that approximation. simpson(F(x),x=0..4,8); simp8:=evalf(%); error8:=actualF-simp8; Pretty good! Spreadsheet summary. We can use the Maple spreadsheet below to see the approximations and errors given by the various methods for approximating the integral of a given function over the interval [NiMlImFH, NiMlImJH]. To use the spreadsheet, enter the function NiMlImZH, the limits of integration NiMlImFH and NiMlImJH, and three values of NiMlIm5H, labeled as NiMlI24xRw==, NiMlI24yRw==, and NiMlI24zRw==, as shown in the next line. f:=sqrt(4+x^4);a:=0;b:=3;n1:=10;n2:=100;n3:=1000; With these values, click the mouse on the blank gray rectangle in the upper left corner of the spreadheet to select all the cells of the spreadsheet and then hit the Enter key. You can see how the spreadsheet was created by clicking on the individual cells and observing the formulas at the top of the page. NiMlIm5H NiMiIzU= NiMiJCsi NiMiJSs1 NiMtJSVsZWZ0RzYjJSJuRw== NiMkIituSz96NSEiKQ== NiMkIis8IT1CPCIhIik= NiMkIis4NC0jPSIhIik= NiMtJSdFX2xlZnRHNiMlIm5H NiMkIipRLCFSNSEiKQ== NiMkIikpUSZ5NSEiKQ== NiMkIigjXCMzIiEiKQ== NiMtJSZyaWdodEc2IyUibkc= NiMkIisrJyp5JkgiISIp NiMkIitdbShSPiIhIik= NiMkIit3bj0lPSIhIik= NiMtJShFX3JpZ2h0RzYjJSJuRw== NiMkISomPidvNyIhIik= NiMkISlYSygzIiEiKQ== NiMkIShyTDMiISIp NiMtJSRtaWRHNiMlIm5H NiMkIitAbyE0PSIhIik= NiMkIitUOTMkPSIhIik= NiMkIisnPS5KPSIhIik= NiMtJSZFX21pZEc2IyUibkc= NiMkIiglZSc+IyEiKQ== NiMkIiZrPiMhIik= NiMkIiQ+IyEiKQ== NiMtJSV0cmFwRzYjJSJuRw== NiMkIitNa1woPSIhIik= NiMkIitNdDkkPSIhIik= NiMkIitYUTUkPSIhIik= NiMtJSdFX3RyYXBHNiMlIm5H NiMkIShISVIlISIp NiMkISZIUiUhIik= NiMkISRTJSEiKQ== NiMtJSVzaW1wRzYjJSJuRw== NiMkIitaRTUkPSIhIik= NiMkIiswTTUkPSIhIik= NiMkIiswTTUkPSIhIik= NiMtJSdFX3NpbXBHNiMlIm5H NiMkIiRlKCEiKQ== NiMkIiIhRiQ= NiMkIiIhRiQ= Let's try another function. f:=x*sin(x);a:=0;b:=2;n1:=10;n2:=100;n3:=1000; NiMlIm5H NiMiIzU= NiMiJCsi NiMiJSs1 NiMtJSVsZWZ0RzYjJSJuRw== NiMkIis5RCoqZjohIio= NiMkIis+eFNCPCEiKg== NiMkIitKRHhSPCEiKg== NiMtJSdFX2xlZnRHNiMlIm5H NiMkIionZSlmIj0hIio= NiMkIikiUSQ9PSEiKg== NiMkIihwJj09ISIq NiMtJSZyaWdodEc2IyUibkc= NiMkIismWzZQIz4hIio= NiMkIis7J3ooZjwhIio= NiMkIis/KDRNdSIhIio= NiMtJShFX3JpZ2h0RzYjJSJuRw== NiMkISomUT9APSEiKg== NiMkISk7Jik9PSEiKg== NiMkISg/Jz09ISIq NiMtJSRtaWRHNiMlIm5H NiMkIituIWY5dSIhIio= NiMkIis7KSplVDwhIio= NiMkIisoMyJmVDwhIio= NiMtJSZFX21pZEc2IyUibkc= NiMkIidMPzghIio= NiMkIiUlRyIhIio= NiMkIiM4ISIq NiMtJSV0cmFwRzYjJSJuRw== NiMkIisqKj4mPXUiISIq NiMkIitvT2ZUPCEiKg== NiMkIitFNmZUPCEiKg== NiMtJSdFX3RyYXBHNiMlIm5H NiMkIScqKjNFISIq NiMkISVvRCEiKg== NiMkISNFISIq NiMtJSVzaW1wRzYjJSJuRw== NiMkIiskNHU6dSIhIio= NiMkIispNCJmVDwhIio= NiMkIisrNmZUPCEiKg== NiMtJSdFX3NpbXBHNiMlIm5H NiMkIiYycSIhIio= NiMkIiIjISIq NiMkIiIhRiQ= Error in Simpson's rule: choosing NiMlIm5H to control error. The error for Simpson's Rule is NiMsJComKSwmJSJiRyIiIiUiYUchIiIiIiZGKComIiQhPUYoKiQpJSJuRyIiJUYoRihGKkYq NiMlImZH ''''(NiMlI211Rw==), where NiMlI211Rw== is some number between NiMlImFH and NiMlImJH. Thus, if |NiMlImZH''''(NiMlInhH)|<=NiMlIk1H for NiMlInhH = NiMlImFH .. NiMlImJH, NiMlIkVH <= NiMqJiksJiUiYkciIiIlImFHISIiIiImRicqJiIkIT1GJyokKSUibkciIiVGJ0YnRik=NiMlIk1H (*) where NiMlIkVH is the absolute value of the error. For NiMvLSUiRkc2IyUieEcqJi0lJXNxcnRHRiYiIiItJSRzaW5HRiZGKw==, used above, we have NiMvJSJhRyIiIQ==, NiMvJSJiRyIiJQ==, and NiMvJSJuRyIiKQ==. We find the fourth derivative of F. a:=0;b:=4;n:=8; F4p:=(D@@4)(F); We look at the graph of this fourth derivative between 0 and 4. plot(F4p(x),x=0..4); By viewing the graph, we could simply choose NiMvJSJNRyIiJg==. But we can be more exact, as follows: M:=evalf(maximize(F4p(x),x=0..4)); Then we compute an upper bound for the absolute error in our approximation. Error:=(b-a)^5/(180*n^4)*M; We note from above that the absolute value of the error was less than this. We can also rewrite the inequality (*) above as NiMlIm5H >= NiMpKiYqJiksJiUiYkciIiIlImFHISIiIiImRiklIk1HRilGKSomIiQhPUYpJSJFR0YpRisqJkYpRikiIiVGKw== so as to choose an appropriate NiMlIm5H so that the absolute error will be less than a chosen value NiMlIkVH. Suppose we wish the error to be less than NiMvJSJFRyQiIiIhIiU=. E:=.0001; n:=evalf(((b-a)^5*M/(180*E))^(1/4)); Thus, since NiMlIm5H must be even and at least as large as this number, choose NiMvJSJuRyIjQw== and use Simpson. simpson(F(x),x=0..4,24); simp24:=evalf(%); error24:=actualF-simp24; We are certainly within the required tolerance.