% function riemann
%
% Computes a Riemann sum and graphs stepwise approximation to function.
% Call is riemann(f,corners,graph) for an inline function f, and
% riemann('f', corners, graph) for a function defined in an mfile.
% corners = [a b c d] is a vector of the corner coordinates of a
% rectangle with corners (a,c), (b,c), (b,d), and (a,d). We assume
% a < b and c < d.
% The third argument is optional. If some number is entered for
% graph, the graph of the stepwise approximation will be displayed.
% If the third argument is omitted from the call, the graph is not
% displayed.
% After the call, user is asked to enter the number n of
% subdivisions in the x direction, and the number m of subdivisions
% in the y direction.
% The Riemann sums are computed using the midpoint of each
% subrectangle.
function out = riemann(f, corners,graph)
a = corners(1); b = corners(2); c = corners(3); d = corners(4);
subdiv = input('enter the number of subdivisions in x and y direction as [n m] ')
n = subdiv(1); m = subdiv(2);
delx = (b-a)/n; dely = (d-c)/m;
x = a:delx:b;
y = c:dely:d;
p = a +.5*delx: delx : b-.5*delx;
q = c +.5*dely: dely : d-.5*dely;
[P,Q] = meshgrid(p,q);
Z = feval(f,P,Q);
disp(' Approximate value of the integral ')
out = sum(sum(Z))*delx*dely;
if nargin < 3
return
end
xplot = zeros(1,2*n);
for j = 1:2:2*n-1
xplot(j) = x((j+1)/2);
end
for j = 2:2:2*n
xplot(j) = a + (j-2)*delx/2 + .99*delx;
end
yplot = zeros(1,2*m);
for i = 1:2:2*m-1
yplot(i) = y((i+1)/2);
end
for i = 2:2:2*m
yplot(i) = c + (i-2)*dely/2 + .99*dely;
end
[Xplot,Yplot] = meshgrid(xplot, yplot);
Zplot = zeros(size(Xplot));
j = 1:2:2*n-1;
k = 2:2:2*n;
for i = 1:m
Zplot(2*i-1,j) = Z(i,:);
Zplot(2*i-1,k) = Z(i,:);
Zplot(2*i, :) = Zplot(2*i-1,:);
end
surf(Xplot, Yplot, Zplot); colormap(cool)