Auxiliary Procedures -- execute but do not modify

> riemann2 := proc(f,xrange,yrange,m,n)

> local A, F, hx, hy, midx, midy, P, rs, x, xlo, xhi, xx, y, ylo, yhi, yy;

> if typematch( xrange,

> 'x'::name='xlo'::realcons..'xhi'::realcons,

> 'matchx' ) then

> assign( matchx );

> hx := (xhi-xlo)/m;

> xx := xlo + [$(0..m)]*hx;

> midx := (xx[1..-2]+xx[2..-1])/2;

> end if;

> if typematch( yrange,

> 'y'::name='ylo'::realcons..'yhi'::realcons,

> 'matchy' ) then

> assign( matchy );

> hy := (yhi-ylo)/n;

> yy := ylo + [$(0..n)]*hy;

> midy := (yy[1..-2]+yy[2..-1])/2;

> end if;

> A := linalg[matrix]( m, n, (i,j)->eval(f, {x=midx[i],y=midy[j]} ) );

> rs := convert( map(op,convert(A,listlist)), `+` ) * hx*hy;

> F := plottools[transform]( (i,j,z)->

> [xlo+(i-1)*hx,ylo+(j-1)*hy,z] );

> P := matrixplot( A, heights=histogram, orientation=[-45,45] );

> return( display( F(P), labels=[x,y,z], title=sprintf("%dx%d Riemann sum = %f", m, n, evalf(rs)) ) );

> end proc:

>

> animate_riemann2 := proc( f, xrange, yrange, N )

> display( [seq(riemann2( f, xrange, yrange, 2^k, 2^k ), k=0..N )],

> insequence=true );

> end proc:

>