%%%%%%%%%% 2020.6.5 %%%%%%%%%
%%%%%%%%%%复化Simpson公式求二重积分%%%%%%%%%
%%%%%%%%%% Liu Deping %%%%%%%%%
clear all;
%%%被积函数及积分上下限导入;
s=input('请输入函数表达式:f = ','s');
f=inline(s);
a = input('请输入积分变量x左边界a的值:');
b = input('请输入积分变量x右边界b的值:');
c = input('请输入积分变量y左边界c的值:');
d = input('请输入积分变量y右边界d的值:');
h1 = input('请输入积分变量x步长h1的值:');
h2 = input('请输入积分变量y步长h2的值:');
m=round((b-a)/h1);
n=round((d-c)/h2);
%%%系数表T,t[i,j]为复化梯形公式的系数;
T=zeros(m+1,n+1);
R=[1,4,1;4,16,4;1,4,1];
for p=0:m/2-1
for q=0:n/2-1
for i=1:3
for j=1:3
T(i+2*p,j+2*q)=T(i+2*p,j+2*q)+R(i,j);
end
end
end
end
fprintf('系数表:');
T
%%%计算各二维节点函数值,并存放于矩阵F中;
F=zeros(m+1,n+1);
for i=0:m
for j=0:n
F(i+1,j+1)=f(a+i*h1,c+j*h2);
end
end
%%%结果输出;
format long;
fprintf('各节点函数值;')
F
fprintf('复化Simpson公式计算结果;')
Snm=(b-a)*(d-c)/(9*m*n)*sum(sum(T.*F))
代码如上,在我们运行代码后发现,命令行窗口中要我们输入函数表达式,输入函数的上下取值范围,以及步长: