在matlab下用梯形法求函数的积分

函数是网上看到的,加上了一点我自己的理解

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

再控制窗口输入一下命令得到的结果为：

Trapezia(0,pi,'sin(x)',eps,7)

***********start*************
k=1 h=0.785398 I=1.570796
k=2 h=0.392699 I=1.896119
k=3 h=0.196350 I=1.974232
k=4 h=0.098175 I=1.993570
k=5 h=0.049087 I=1.998393
k=6 h=0.024544 I=1.999598
k=7 h=0.012272 I=1.999900

***********end*************

可以看出这个函数能都根据要就精确度E来计算结果，但是函数不断的重复运算使运算时间比较长。