> 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);