% FILE hlp052.m begins.
function hlp052
a=0; b=4; c=0;d=1;
fprintf('Using doubleint0 (Riemann sum):\n')
exact=  (d^2- c^2)*(b^2 -a^2)/4;
for m=20:40:100
for n=20:40:100
        numer=doubleint0(a,b,c,d,m,n);
        fprintf('  %12.3e', numer-exact);
end
fprintf('\n')
end

% REMARK: The built-in MATLAB function gives
% accuracy of the order 10E-16
% r2=dblquad('x.*y',a,b,c,d);
% exact=  (d^2- c^2)*(b^2 -a^2)/4;
% fprintf('%12.3e\n', exact-r2)
exact=  (d^2- c^2)*(b^2 -a^2)/4;
fprintf('Using doubleint2 (Simpson Rule):\n')
for m=20:40:100
for n=20:40:100
        numer=doubleint2(a,b,c,d,m,n);
        fprintf('  %12.3e', numer-exact);
end
fprintf('\n')
end

function z= myf(x,y)
z= x.*y;



function val=doubleint0(x1,x2,y1,y2,m,n)
% Riemann Sum
a=x1; b=x2;
s=0.0;
dx=(b-a)/m;
dy=(y2-y1)/n;
for ii=0:(m-1)
     x=a+ii*dx;
for jj=0:(n-1)
     y=y1+jj*dy;
     s=s+myf(x,y);  
end
end
val=s*dx*dy;

function w=doubleint2(x1,x2,y1,y2,m,n)
% Note calling program must supply myf  
% Compute \int_x1^x2  \int_{y1}^{y2} myf(x,y) dx dy
% This is based on the one-dimensional Simpson rule
% \int_a^b f(x) dx = (h/3) [y_0+ 4y_1 + 2y_2+ 4y_3 + ...+2y_{2n-2}+ 
%                      4 y_{2n-1}+ y_{2n}] 
% x_k= a+ k*(b-a)/(2n),   y_k=f(x_k), k=0,..., 2n

a=x1; b=x2;
s=0.0;
wx=ones(1,1+m);
wx(2:2:(end-1))=4*ones(size(wx(2:2:(end-1))));
wx(3:2:(end-2))=2*ones(size(wx(3:2:(end-2))));
wx=(1/3)*wx;
wy=ones(1,1+n);
wy(2:2:(end-1))=4*ones(size(wy(2:2:(end-1))));
wy(3:2:(end-2))=2*ones(size(wy(3:2:(end-2))));
wy=(1/3)*wy;

for ii=0:(m)
     x=a+ii*(b-a)/m;
     dy=(y2-y1)/n;
for jj=0:(n)
     y=y1+jj*dy;
     s=s+wx(ii+1)*wy(jj+1)*myf(x,y)*dy*(b-a)/m;
end
end
w= s;
% FILE hlp052.m ends.

