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:
>