Built-In Mathematica Command to Calculate Integrals

In Sections 5.3-5.5 we will learn how to calculate certain types of integrals by hand, without resorting to interpreting them as areas under curves. We will want to be able to check our answers, and we can do that using Mathematica. For example, consider the following integral we want to calculate:

(∫_(-1)^3)(x^2)dx

Mathematica can do this a number of ways:

Integrate[x^2,{x,-1,3}]

(∫_(-1)^3)x^2\[DifferentialD]x

NIntegrate[x^2,{x,-1,3}]

Note that the varible x truly is a dummy index, since

(∫_(-1)^3)y^2\[DifferentialD]y
(∫_(-1)^3)w^2\[DifferentialD]w
(∫_(-1)^3)zipper^2\[DifferentialD]zipper
(∫_(-1)^3)θ^2\[DifferentialD]θ

Using Mathematica  to Plot Those Rectangles in Riemann Sum

(∫_a^b)f(x)dx=(lim_(n→∞))(∑_(i=1)^n)f((SuperStar[(Subscript[x, i])]))(((b-a))/n)~(∑_(i=1)^n)f((SuperStar[(Subscript[x, i])]))(((b-a))/n) if n is large.

Using Mathematica to plot the rectangles which we use to estimate the area under a curve requires a small amount of work to set the plot up. However, this work is well worth it! Remember, we worked out in class the following:

Right Hand Side Endpoints in the partition: (Subscript[R, n])

We are given a function that we wish to find the area under, using a finite number of rectangles. The xlimits of the integration a and b are given, as well as the number of intervals (partitions) to use. For example, we might be given:

a=1;
b=3;
n=10;
f[x_] = 1/Sqrt[x^3+1]

We can plot the function:

plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}]

Here is how we can generate the rectangles.

deltax = (b-a)/n;
data = Table[Graphics[{Orange,Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+i*deltax]}]}],{i,1,n}];

This creates a data set with images of the rectangles. The command in red is what creates the rectangle. The points are the {(Subscript[x, min]),(Subscript[y, min])} and {(Subscript[x, max]),(Subscript[y, max])}, which for using right hand endpoints as the height of the rectangle means {a+(i-1)*deltax,0},{a+i*deltax,f[a+i*deltax]}. The Blue, Opacity, and EdgeForm are just style options to make the rectangles look better.

Let's take a look at the data set:

data

The scaling of the rectangles is not apparent above, but if we plot them all together we see that they are the correct scale.

Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

Let's put this all together so that it easier to work with:

a=1;
b=3;
n=10;
f[x_]=(1/(Sqrt[(x^3)+1]));
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);data = Table[Graphics[{Orange,Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+i*deltax]}]}],{i,1,n}];
Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

We can easily change the number of partitions, the function, or anything else we desire:

a=-1;
b=2;
n=25;
f[x_]=(1/(Sqrt[(x^3)+1]));
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);data = Table[Graphics[{Orange,Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+i*deltax]}]}],{i,1,n}];
Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

Here is how we can add up all the areas in the rectangles, and approximate the integral:

R[n_,a_,b_] := N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*i]]

TableForm[Table[{n,R[n,-1,2]},{n,1,25}]]

Net Area: Colour the rectangles Blue and Orange

Our textbook colours the rectangles blue when they are below the x-axis, and orange when they are above the x-axis. This is done since we subtract the area below the axis from the area above to give us the net area. We can do that by defining a simple function that chooses the correct color for us!

colour[x_]:=Piecewise[{{Blue,x<0},{Orange,x≥0}}]

Here is the sketch using Right endpoints, with our new colour selector:

a=-1;
b=3;
n=12;
f[x_]=Cos[x];
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);data = Table[Graphics[{colour[f[a+i*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+i*deltax]}]}],{i,1,n}];
Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

Left Hand Side Endpoints in the partition: (Subscript[L, n])

We can easily change the code to work with left hand endpoints. All we have to do is change the maximum y value in the Rectangle command to be the left endpoint, not the right.

Here is the sketch using Right endpoints:

colour[x_]:=Piecewise[{{Blue,x<0},{Orange,x≥0}}]
a=-1;
b=3;
n=12;
f[x_]=Cos[x];
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);data = Table[Graphics[{colour[f[a+i*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+i*deltax]}]}],{i,1,n}];
Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

And here it is using the Left endpoints (the only change is in red below, where we evaluate at the left endpoint of the partition):

colour[x_]:=Piecewise[{{Blue,x<0},{Orange,x≥0}}]
a=-1;
b=3;
n=12;
f[x_]=Cos[x];
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);data = Table[Graphics[{colour[f[a+(i-1)*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+(i-1)*deltax]}]}],{i,1,n}];
Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

Here is how we can add up all the areas in the rectangles, and approximate the integral using left endpoints:

L[n_,a_,b_] := N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*(i-1)]]

TableForm[Table[{n,L[n,-1,2]},{n,1,20}]]

Midpoints in the partition: (Subscript[M, n])

A similar modification allows us to use the midpoint of each partition

colour[x_]:=Piecewise[{{Blue,x<0},{Orange,x≥0}}]
a=-1;
b=3;
n=12;
f[x_]=Cos[x];
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);data = Table[Graphics[{colour[f[a+((b-a)/n)*(i-(1/2))]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+((b-a)/n)*(i-(1/2))]}]}],{i,1,n}];
Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

Here is how we can add up all the areas in the rectangles, and approximate the integral:

M[n_,a_,b_] := N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*(i-(1/2))]]

TableForm[Table[{n,M[n,-1,2]},{n,1,20}]]

Animation: (Subscript[L, n]) as we vary n

colour[x_]:=Piecewise[{{Blue,x<0},{Orange,x≥0}}]

By simply wrapping a Manipulate command around the plot, we can generate animations of  how (Subscript[L, n])changes as the value of n increases.

Manipulate[
a=-1;
b=3;
f[x_]=Cos[x];
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);data = Table[Graphics[{colour[f[a+(i-1)*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+(i-1)*deltax]}]}],{i,1,n}];
Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]
,{n,10,50,1}]

Here is another example.

Manipulate[
a=0;
b=3;
f[x_]= x^3-6x;
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);data = Table[Graphics[{colour[f[a+(i-1)*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+(i-1)*deltax]}]}],{i,1,n}];
Show[plot1,data,ImageSize->{400,Automatic},AxesLabel->{"x","y"}]
,{n,10,50,1}]

A More Complicated Animation, including Riemann Sums

Here we plot (Subscript[L, n])and (Subscript[R, n]) side by side, and include the value of the Riemann sums.

Here is how we can calculate Riemann sums in Mathematica:

R[n_,a_,b_] :=  N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*i]]
L[n_,a_,b_] := N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*(i-1)]]

Manipulate[
a=0;
b=3;
f[x_]= x^3-6x;
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);dataL = Table[Graphics[{colour[f[a+(i-1)*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+(i-1)*deltax]}]}],{i,1,n}];
dataR = Table[Graphics[{colour[f[a+i*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+i*deltax]}]}],{i,1,n}];
plotR=Show[plot1,dataR,ImageSize->{400,Automatic},AxesLabel->{"x","y"}];
plotL=Show[plot1,dataL,ImageSize->{400,Automatic},AxesLabel->{"x","y"}];
Show[GraphicsGrid[{{plotL,plotR},{L[n,a,b],R[n,a,b]}}],ImageSize->{500}]
,{n,10,50,1}]

A Difficult Integral

Here is an example of a difficult integral. We haven't looked at evaluating integrals yet, but let's just see if Mathematica can do it!

(∫_(-2)^3)Exp[-(x^4+x^2)/10]\[DifferentialD]x

Mathematica could always do every derivative we asked it to, no matter how complicated. The fact the Mathematica could not do this integral is the first indication that integration is a more complicated process than differentiation.

Although we could not get a numerical value for this integral, we can still use our Riemann sums to get an approximation to the integral.

Let's see what our animation can tell us:

R[n_,a_,b_] :=  N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*i]]
L[n_,a_,b_] := N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*(i-1)]]

Manipulate[
a=-2;
b=3;
f[x_]=Exp[-(1/10) ((x^4)+(x^2))];
plot1 = Plot[f[x],{x,a,b},PlotStyle -> {Blue,Thick},AxesOrigin->{0,0}];
deltax=((b-a)/n);dataL = Table[Graphics[{colour[f[a+(i-1)*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+(i-1)*deltax]}]}],{i,1,n}];
dataR = Table[Graphics[{colour[f[a+i*deltax]],Opacity[0.2],EdgeForm[Thick],Rectangle[{a+(i-1)*deltax,0},{a+i*deltax,f[a+i*deltax]}]}],{i,1,n}];
plotR=Show[plot1,dataR,ImageSize->{400,Automatic},AxesLabel->{"x","y"}];
plotL=Show[plot1,dataL,ImageSize->{400,Automatic},AxesLabel->{"x","y"}];
Show[GraphicsGrid[{{plotL,plotR},{L[n,a,b],R[n,a,b]}}],ImageSize->{500}]
,{n,2,100,1}]

It looks like the area under the curve is about 2.88.

Mathematica can verify this for us if we use the NIntegrate command

NIntegrate[Exp[-(x^4+x^2)/10],{x,-2,3}]

A Closer Look for Very Large Values of  n

Like with derivatives, numerical round-off errors can cause problems with these numerical techniques to do intgerals.

Here is a purely numeric example:

a=0.0;
b=3.0;
f[x_]=Cos[x]Sin[x]Exp[x];
R[n_,a_,b_] :=  N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*i]]
L[n_,a_,b_] := N[(((b-a))/n) (∑_(i=1)^n)f[a+((b-a)/n)*(i-1)]]

TableForm[Table[{n,L[n,a,b],R[n,a,b]},{n,1,100}]]

We see that we seem to be approaching -4.2, which is close to what Mathematica tells us the integral is:

Integrate[f[x],{x,a,b}]

Here is a plot that shows that (Subscript[L, n]) is approaching -4.2.

list=Table[{n,L[n,a,b]},{n,1000,10000,1000}];
ListPlot[list,Joined->False,ImageSize->{500,Automatic},AxesLabel->{"n","(Subscript[L, n])"},BaseStyle->Large,PlotStyle->PointSize[0.02]]

So we might expect that choosing super large values for n would give us numbers even closer to -4.2:

TableForm[Table[L[10^n,a,b],{n,1,5}]]

This was a bit slow, but it would be even slower except for the fact that Mathematica is pretty good at manipulating summations. Look at what Mathematica is using for our L[m,a,b]:

L[m,a,b]

That's not what we typed in! Mathematica has simplified the summation into something else, that is faster to evaluate and suffers from less numerical imprecision than the original form.

This ability to simplify the summation is what you would have to do to be able to evaluate integrals from the definition, which is one of the reasons why the Fundamental Theorem of Calculus is so important--by connecting integrals and derivatives, we avoid having to do this. Yay!