volumesim :=
proc(f,a,b,c,d,m,numpts)
#probabilistic simulation of volume under z = f(x,y) from x=a to x=b and y=c and y=d:
# NOTE: z = f(x,y) must be continuous and >= 0 on x element of [a,b] and y element of [c,d]:
#input params:
# f = the z = f(x,y) function, defined in Maple worksheet:
# a,b = x limits of integration:
# c,d = y limits of integration:
# m = 2: the maximum value of f(x,y) for x in its bounds and y in its bounds:
# numpts = number of random points generated for simulation:
#output: estimate of volume under curve:
local r,ct,ctunder,xi,yi,zi,approxvolume;
r := rand(1..1000); #random number betw 1 & 1000:
ctunder := 0;
for ct from 1 to numpts #generate random pts in box &:
do #count how many are under surface:
xi := evalf( (b-a) * (r()/1000) + a );
yi := evalf( (d-c) * (r()/1000) + c );
zi := evalf(m * (r()/1000) );
if zi <= f(xi,yi) then ctunder := ctunder + 1 fi;
od;
approxvolume := m * (b-a) * (d-c) * (ctunder/numpts); #estimate the volume:
print(evalf(approxvolume));
end;