% Experiment 17.2: randomly chosen Riemann sums.
%
% The program computes Riemann sums for the function f(x,y) = x^2+y^2 on
% the rectangle [0,1]x[0,1], plots the grid G with intermediate points
% and the approximation of the graph by the top faces of the cuboids.
%
% Number of summands: n; the grid is generated by equidistant partitions
% Z_x = {0, deltax, 2*deltax, ..., 1}, Z_y = {0, deltay, 2*deltay, ..., 1} 
% with deltax = deltay = 1/n.
%
% The intermediate points are chosen with the help of 2*n^2 random numbers,
% uniformly distributed in each rectangle of the grid.
%
% Call of program after choice of n by 
%
%       mat17_1(n);
%
% (n <= 50 is recommended).

function fct(n);

% defaults
if nargin < 1 
   n = 10;
end
if n > 50
  'To limit computing time, n <= 50 is recommended - or else comment Figure(2)'
else
   
% partitions
deltax = 1/n;
deltay = 1/n;
Zx = 0:deltax:1-deltax;
Zy = 0:deltay:1-deltay;

% random intermediate points and corresponding function values
Rx = rand(n,n);
Ry = rand(n,n);

xi = zeros(size(Zx));
xi = zeros(size(Zy));
f = zeros(n,n);
for i = 1:n
    for j=1:n
        xi(i,j) = Zx(i) + deltax*Rx(i,j);
        eta(i,j) = Zy(j) + deltay*Ry(i,j);
        f(i,j) = xi(i,j)^2 + eta(i,j)^2;
    end
end

% value of Riemann sum
I = sum(sum(f))*deltax*deltay

% plot of grid and intermediate points
figure
hold on
Zx1 = [Zx 1];
Zy1 = [Zy 1];
for k = 1:n+1
        plot([Zx1(k) Zx1(k)],[0 1],'k')
        plot([0 1],[Zy1(k) Zy1(k)],'k')
end
plot(xi,eta,'r.')
title('Grid and intermediate points')
hold off


% 3D-image of top faces of the cuboids
figure
hold on
for i = 1:n
    for j = 1:n
 % top face
 fill3([Zx1(i) Zx1(i+1) Zx1(i+1) Zx1(i)], [Zy1(j) Zy1(j) Zy1(j+1) Zy1(j+1)], [f(i,j) f(i,j) f(i,j) f(i,j)],[1 1 1 1])
 % bottom grid
 fill3([Zx1(i) Zx1(i+1) Zx1(i+1) Zx1(i)], [Zy1(j) Zy1(j) Zy1(j+1) Zy1(j+1)], [0 0 0 0],'w')
    end
end

% boundary of the surface z = x^2 + y^2
t = 0:0.01:1;
plot3(t, zeros(size(t)), t.^2, 'k', 'Linewidth',1.5)
plot3(zeros(size(t)), t, t.^2, 'k', 'Linewidth',1.5)
plot3(t, ones(size(t)), 1 + t.^2, 'k', 'Linewidth',1.5)
plot3(ones(size(t)), t, 1 + t.^2, 'k', 'Linewidth',1.5)
plot3([1 1],[0 0],[0 1],'k','Linewidth',1)
plot3([0 0],[1 1],[0 1],'k','Linewidth',1)
plot3([1 1],[1 1],[0 2],'k','Linewidth',1)
title('Top faces of cuboids')
view([-40 60])

hold off
end