erchongjifen.m

该程序为重积分的计算。应用复化梯形公式

function ex5_23_1
syms x y f;
f=input('请输入被积函数f(x,y)：');
a=input('请输入X的下限：');
b=input('请输入X的上限：');
c=input('请输入Y的下限：');
d=input('请输入Y的上限：');
e=input('请输入精度：');

h=b-a;
k=d-c;
T0=zeros(1,5);
T1=zeros(1,4);
T2=zeros(1,3);
T3=zeros(1,2);

IT=0;%积分结果

F1=fun(f,a,c);
F2=fun(f,a,d);
F3=fun(f,b,c);
F4=fun(f,b,d);
T0(1)=h*k*(F1+F2+F3+F4)/4;

p=1; %二分次数
while p<=4
    h=h/2;
    k=k/2;
    m=(b-a)/h;%x区间等分数
    n=(d-c)/k;%y区间等分数

    
    Xx=zeros(1,m+1);
    Yy=zeros(1,n+1);
    
    for i=0:m
        Xx(i+1)=a+i*h;
    end
    for j=0:n
        Yy(j+1)=c+j*k;
    end
    
    for i=1:m      
        for j=1:n
            F1=fun(f,Xx(i),Yy(j));
            F2=fun(f,Xx(i),Yy(j+1));
            F3=fun(f,Xx(i+1),Yy(j));
            F4=fun(f,Xx(i+1),Yy(j+1));
            T0(p+1)=T0(p+1)+h*k*(F1+F2+F3+F4)/4;
        end
    end
    
    p=p+1;
end
 

for q=1:4
    T1(q)=(4*T0(q+1)-T0(q))/3;
end
if abs(T1(2)-T1(1))<e
    IT=T1(2);   
else for q=1:3
        T2(q)=(16*T1(q+1)-T1(q))/15;
    end
    if abs(T2(2)-T2(1))<e
         IT=T2(2);
    else for q=1:2
            T3(q)=(64*T2(q+1)-T2(q))/63;
        end
        if abs(T3(2)-T3(1))<e
           IT=T3(2);
        else IT=(256*T3(2)-T3(1))/255;
        end
    end
end

fprintf('\n积分结果是：%f',IT)
     

function F=fun(f,m,n)
syms x y;
F=subs(f,{x,y},{m,n});