trapezia.m

用于求积分

%用来就数值积分
%fx是由syms定义的函数表达式
function Trapezia(a,b,fx,E,N)
fprintf('\n***********start*************\n');
h=(b-a)/2;
a=a+E/10;%用来减少出现除数为零的情况
               %比如对含有1/x的函数积分且积分区间从0开始
               %在梯形公式及Simpson公式算法中必然出现除数为零的情况
               %当然我们可以在这之前就对被积函数做一定变形处理 比如分部积分将分母移走
               %不过我更愿意采用更简单的做法 将积分下限加一个非常小的数 使之不等于零
               %这种处理方法对于精度要求几乎没有影响
               %当然who有更好的方法欢迎交流一下！
T0=h*(feval(fx,a)+feval(fx,b));
T=0;
for m=1:N
    F=0;
    for k=1:2^(m-1)
        F=F+feval(fx,a+(2*k-1)*h);
    end
    T=0.5*T0+h*F;
    if abs(det(T-T0))<3*E
        break;
    else
        h=h/2;
        T0=T;
    end
     fprintf('k=%d\th=%f\tI=%f\n',m,h,T);
end
fprintf('\n***********end*************\n');
end