数值分析上机题 5 2010/12/6 一、题目 P.236, 23.（上机题）重积分的计算 二、程序 38 - (1) 计算 I(f)的通用程序（Matlab）： 表 1 % Get the approximation of a double integral of f(x,y) % on a rectangular region R={(x,y)|a<=x<=b,c<=y<=d} % with given step m=(b-a)/h and n=(d-c)/k % and given error epsilon function [bestres,res] = getDoubleIntegral(f,a,b,c,d,m,n,epsilon) res = []; count = 1; res(count,1) = getT(f,a,b,c,d,2^count*m,2^count*n); count = 2; res(count,1) = getT(f,a,b,c,d,2^count*m,2^count*n); res(count-1,2) = 4/3*res(count,1) - 1/3*res(count-1,1); res1 = res(1,1); res2 = res(2,1); co = [4/3 -1/3;16/15 -1/15;64/63 -1/63]; while abs(res1-res2) > epsilon count = count + 1; res(count,1) = getT(f,a,b,c,d,2^count*m,2^count*n); %#ok
for i = 1:size(co,1) if count > i res(count-i,i+1) = co(i,1)*res(count-i+1,i) + co(i,2)*res(count-i,i); %#ok end end i = size(co,1); while size(res,1) - i - 1 < 0 i = i - 1; end res1 = res(count-i,i); res2 = res(count-i+1,i); bestres = res2; end end % Get approximation of a double integral of f(x,y) % on a rectangular region R={(x,y)|a<=x<=b,c<=y<=d} % with given steps m=(b-a)/h and n=(d-c)/k % using The Composite Trapezoid Rule function res = getT(f,a,b,c,d,m,n) h = (b-a)/m; k = (d-c)/n; x = a + (0:m)'*h; y = c + (0:n)'*k; res = 0; for i = 1:m for j = 1:n res = res + h*k/4 * (f(x(i),y(j)) + f(x(i),y(j+1)) + f(x(i+1),y(j)) + f(x(i+1),y(j+1))); end end end 3 8 - (2 ) 求解： 表 2 function exp5 f = @myfun; a = 0; b = pi/6; c = 0; d = pi/3; m = 1; n = 1; epsilon = 0.5e-5; [bestres,res] = getDoubleIntegral(f,a,b,c,d,m,n,epsilon); fprintf('================================================================================\n') fprintf('Result is %0.7f, detail:\n',bestres) fprintf('======...
