Page 253 - Applied Numerical Methods Using MATLAB

242 NUMERICAL DIFFERENTIATION/ INTEGRATION
(cf) The MATLAB built-in routine dblquad() can accept the boundaries of integration
region only given as numbers. Therefore, if we want to use the routine in comput-
ing a double integral for a nonrectangular region D, we should define the integrand
function f (x, y) for a rectangular region R ⊇ D (containing the actual integration
region D) insuch a waythat f (x, y) = 0for (x, y) /∈ D; that is, the value of the
function becomes zero outside the integration region D, which may result in more
computations.

function INTfxy = int2s(f,a,b,c,d,M,N)
%double integral of f(x,y) over R = {(x,y)|a <= x <= b, c(x) <= y <= d(x)}
% using Simpson’s rule
if ceil(M) ~= floor(M) %fixed width of segments on x
hx = M; M = ceil((b - a)/hx);
end
if mod(M,2) ~= 0, M=M+1;end
hx = (b - a)/M; m = 1:M+1;x=a+(m- 1)*hx;
if isnumeric(c), cx(m) = c; %if c is given as a constant number
else cx(m) = feval(c,x(m)); %in case c is given as a function of x
end
if isnumeric(d), dx(m) = d; %if c is given as a constant number
else dx(m) = feval(d,x(m)); %in case d is given as a function of x
end
if ceil(N) ~= floor(N) %fixed width of segments on y
hy = N; Nx(m) = ceil((dx(m)- cx(m))/hy);
ind = find(mod(Nx(m),2) ~= 0); Nx(ind) = Nx(ind) + 1;
else %fixed number of subintervals
if mod(N,2) ~= 0, N = N +1; end
Nx(m) = N;
end
form=1:M+1
sx(m) = smpsns fxy(f,x(m),cx(m),dx(m),Nx(m));
end
kodd = 2:2:M; keven = 3:2:M - 1; %the set of odd/even indices
INTfxy = hx/3*(sx(1) + sx(M + 1) + 4*sum(sx(kodd)) + 2*sum(sx(keven)));
function INTf = smpsns fxy(f, x, c, d, N)
%1-dimensional integration of f(x,y) for Ry = {c <= y <= d}
if nargin < 5, N = 100; end
if abs(d - c)< eps | N <= 0, INTf = 0; return; end
if mod(N,2) ~= 0, N=N+1;end
h = (d - c)/N; y = c+[0:N]*h; fxy = feval(f,x,y);
fxy(find(fxy == inf)) = realmax; fxy(find(fxy == -inf)) = -realmax;
kodd = 2:2:N; keven = 3:2:N - 1; %the set of odd/even indices
INTf = h/3*(fxy(1) + fxy(N + 1) + 4*sum(fxy(kodd)) + 2*sum(fxy(keven)));
%nm510: the volume of a sphere
x = [-1:0.05:1]; y = [0:0.05:1]; [X,Y] = meshgrid(x,y);
f510 = inline(’sqrt(max(1 - x.*x - y.*y,0))’,’x’,’y’);
Z = f510(X,Y); mesh(x,y,Z);
a=-1;b=1;c =0; d=inline(’sqrt(max(1 - x.*x,0))’,’x’);
Vs1 = int2s(f510,a,b,c,d,100,100) %with fixed number of segments
error1 = Vs1 - pi/3
Vs2 = int2s(f510,a,b,c,d,0.01,0.01) %with fixed segment width
error2 = Vs2 - pi/3