m5_4.asv

matlab 相关的数值积分 例如高斯积分 龙贝格积分等等 可参考。

%%%%%%%%  Double Integration by the Simpson's Rule  %%%%%%%%
clear;
nx=input('Number of intervals in x_direction:');
ny=input('Number of intervals in y_direction:');
a=input('Lower limit of x,a?');
b=input('Lower limit of x,b?');
hx=(b-a)/nx;
x=[a:hx:b];
cy=c(x);
dy=d(x);
hy=(d(x)-c(x))/ny;
Result=0;
temp=0;
yy=zeros(1,ny);

%%%%%%%%  G(x_i)  %%%%%%%%%%%%%%%%%%
for i=1:nx+1
         y=[cy(i):hy(i):dy(i)];
         for j=1:ny
            yy(j)=y(j)+hy(i)/2;  %%% 中点 %%%%%%
         end
         zz=f(x(i),yy);
         z=f(x(i),y);
         temp=hy(i)/6*(4*sum(zz)+2*sum(z)-f(x(i),cy(i))-f(x(i),dy(i)));
         if i>1&i<nx+1 temp=2*temp;end
         fprintf('%12')
         Result=Result+temp;
end

%%%%%%%%  G(x_i-1/2)  %%%%%%%%%%%%%%%%%%
yy=zeros(1,ny);
for i=1:nx
    xx(i)=x(i)+hx/2;
end

for i=1:nx
         y=[cy(i):hy(i):dy(i)];
         for j=1:ny
            yy(j)=y(j)+hy(i)/2;
         end
         zz=f(xx(i),yy);
         z=f(xx(i),y);
         temp=hy(i)/6*(4*sum(zz)+2*sum(z)-f(xx(i),cy(i))-f(xx(i),dy(i)));
         temp=4*temp;
         Result=Result+temp;
end

Result=hx/6*Result