>
Integrals and Riemann Sums in One, Two, and Three Variables
Approximating sums in dimension one
Here are some examples of commands that use Maple to calculate approximating
sums for integrals in one, two, and three variables.
The first few inputs may look mysterious; they simply define some of the commands
we'll need for future reference. There's no need to worry too much about
the precise details of how they're defined, but do watch the syntax in the
examples.
> midpointsum := proc (f,interval,n) local a,b,x,i,step,total;
> a := op(1,op(2,interval));
> b := op(2,op(2,interval));
> x := op(1,interval);
> step := (b-a)/n;
> total := 0;
> for i from 1 to n do
> total := total + evalf( subs(x=a+i*step-step/2,f) );
> od;
> total := evalf( total*step );
> end:
>
The previous command defined a midpoint rule Riemann sum for functions of
one variable. Here are some examples to show how it works. The
last input tells how many subdivisions of the interval of integration
the rule should use. We'd expect that with more subdivisions, the answer ought
to get closer and closer to the `exact' value of the integral.
Suppose, e.g., that we want to integrate sin(x) from x=0 to x=1. The answer
is easy to find by antidifferentiation:
> int( sin(x),x=0..1);
Let's find a numerical version:
> evalf(%);
Now let's compare this with some approximations to the same integral,
using the midpoint rule. First we'll use just four subdivisions:
> midpointsum( sin(x), x=0..1, 4);
The answer isn't too bad an approximation, given the small number of subdivisions.
But let's bump that number up.
> midpointsum( sin(x), x=0..1, 8);
> midpointsum( sin(x), x=0..1, 100);
The estimates seem to behave as they should---larger numbers of subdivisions
tend to give better estimates. But it's striking how well the method works even
for small numbers of subdivisions.
Try the same idea on some other functions and/or other intervals.
Approximating sums in dimension two
Now let's do some approximating sums for a function f(x,y) of two variables. The
idea is to approximate the integral of such a function over a rectangle [a,b] x [c,d] .
First we have to tell Maple the two-dimensional version of the midpoint sum
rule. (Think about what midpoint means in this case!) As above, there's no
need to worry about details of the Maple procedure below---what it does
is what's important.
> doublemidpointsum := proc (f,interval1,interval2,n)
> local a,b,c,d,x,y,i,j,stepx,stepy,total;
> a := op(1,op(2,interval1));
> b := op(2,op(2,interval1));
> c := op(1,op(2,interval2));
> d := op(2,op(2,interval2));
> x := op(1,interval1);
> y := op(1,interval2);
> stepx := (b-a)/n;
> stepy := (d-c)/n;
> total := 0;
> for i from 1 to n do
> for j from 1 to n do
> total := total +
> evalf( subs(x=a+i*stepx-stepx/2,y=c+j*stepy-stepy/2,f));
> od;
> od;
> total := evalf( total*stepx*stepy );
> end:
Notice that there are now two intervals involved---one in the x-variable and
one in the y-variable.
Let's apply the rule to the function f(x,y) = sin(x)sin(y), on the rectangle
[0,1] x [0,1], using 10 subdivisions in each direction. (This produces a grid
of 100 squares.)
> doublemidpointsum( sin(x)*sin(y), x=0..1, y=0..1, 10);
The answer is an approximation to the integral of sin(x)sin(y) over the same
square. Geometrically, such an integral tells the volume of the solid whose base
is on the square [0,1] x [0,1] in the xy-plane, and whose ``top'' is the surface
z=f(x,y). Let's look at this figure, to see whether the number above seems
like any sort of sensible approximation to its volume:
> plot3d( sin(x)*sin(y), x=0..1, y=0..1, axes=boxed);
Take a close look. Does it seem possible that the volume is around 0.2?
To find the volume ``exactly,'' we can use an iterated integral, i.e., we can
integrate in one variable and then the other:
> int( sin(x)*sin(y), x=0..1);
The result is a function of y (as it should be ---think about it).
Now we can integrate in the other variable:
> int(%,y=0..1);
We could have done the same thing all in one step, if we'd wanted:
> int( int( sin(x)*sin(y), x=0..1), y=0..1 );
Let's get a decimal value:
> evalf(%);
Notice how similar that number is to the one we found using the midpoint rule. With
fewer subdivisions, we get less good agreement (but still not too bad):
> doublemidpointsum( sin(x)*sin(y), x=0..1, y=0..1, 4);
Try the same idea on some other functions f(x,y) and some other intervals.
Approximating sums in dimension three
Approximating sums work much the same way for a function g(x,y,z) of three variables.
This time, we approximate the integral of such a function over a cube
[a,b] x [c,d] x [e,f] .
Here is Maple's three-dimensional version of the midpoint sum. It's even
messier than the earlier versions, and you can safely ignore its Maple code:
> triplemidpointsum := proc (g,interval1,interval2,interval3,n)
> local a,b,c,d,e,f,x,y,z,i,j,k,stepx,stepy,stepz,total;
> a := op(1,op(2,interval1));
> b := op(2,op(2,interval1));
> c := op(1,op(2,interval2));
> d := op(2,op(2,interval2));
> e := op(1,op(2,interval3));
> f := op(2,op(2,interval3));
> x := op(1,interval1);
> y := op(1,interval2);
> z := op(1,interval3);
> stepx := (b-a)/n;
> stepy := (d-c)/n;
> stepz := (f-e)/n;
> total := 0;
> for i from 1 to n do
> for j from 1 to n do
> for k from 1 to n do
> total := total + evalf( subs(x=a+i*stepx-stepx/2,
> y=c+j*stepy-stepy/2,
> z=e+k*stepz-stepz/2,g));
> od; od; od;
> total := evalf( total*stepx*stepy*stepz );
> end:
>
Let's use the command to approximate the integral of the constant function g(x,y,z)=x+y+z
over the cube [0,1] x [0,2] x [0,3] . Note that setting n=4 here really amounts
to using 4^3 = 64 subdivisions. (Setting n=10 would involve 1000 subdivisions!)
> triplemidpointsum( x+y+z, x=0..1, y=0..2, z=0..3, 4);
That number can't readily be interpreted geometrically, but it does agree with what
we'd get by iterated integration:
> int( int( int( x+y+z,x=0..1), y=0..2), z=0..3);