执行函数为 mymulNewtonCotes.m 1、使用方法： Step1：在 MATLAB 命令窗口输入被积函数 输入应为：ft=@(t)t.*exp(t^2/2)。 2 t 2 t e dt 。 1  0 Step2：执行函数。输入形式为 mymulNewtonCotes(ft,a,b,m,n)； 其中ft—被积函数，此体重ft=@(t)t.*exp(t^2/2)，已经在第一步赋值； a—积分下限，本题中为0； b—积分上限，本题中为1； m—将区间[a,b]等分的子区间数量，本题可选为10； n—采用的Newton-Cotes公式的阶数，必须满足n<8，否则积分没法 保证稳定性。 当n=1时，即为复化梯形公式；n=2时，即为复化复化辛普森公式。 所以，分别输入mymulNewtonCotes(ft,0,1,10,1)和 mymulNewtonCotes(ft,0,1,10,2)就可以得到两种方法的积分计算结果。 2、计算结果 而根据积分运算，可得： 2 t 2 t e dt  1  0 1  0 2 t e d 2 ( )  1 2 0  2 t 2 x e dx x  e 1 2  e 1 2 0 0  e  1.6487 1 0.6487   说明复化梯形和复化辛普森公式计算出的结果基本一致，与实际结果相符。 

3、程序代码 function yy = mymulNewtonCotes(ft,a,b,m,n) % 复化Newton-Cotes数值积分公式，即在每个子区间上使用Newton-Cotes公式，然后求和, % 参考的输入形式为mymulNewtonCotes(ft,0,1,10,2) % 参数说明: % ft——被积函数，此题中ft=@(t)t.*exp(t^2/2) % a——积分下限 % b——积分上限 % m——将区间[a,b]等分的子区间数量 % n——采用的Newton-Cotes公式的阶数，必须满足n<8，否则积分没法保证稳定性 % % (1)n=1时为复化梯形公式 (2)n=2时为复化辛普森公式 xx = linspace(a,b,m+1); for l = 1:m s(l) = myNewtonCotes(ft,xx(l),xx(l+1),n); end yy = sum(s); function [y,Ck,Ak] = myNewtonCotes(ft,a,b,n) % 牛顿-科特斯数值积分公式 % Ck——科特斯系数 % Ak——求积系数 % y——牛顿-科特斯数值积分结果 xk = linspace(a,b,n+1); for j = 1:n+1 ff(j) = ft(xk(j)); end % 计算科特斯系数 for i=1:n+1 

k=i-1; Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k),0,n); end % 计算求积系数 Ak=(b-a)*Ck; % 求和算积分 y=Ak*ff'; function f=intfun(t,n,k) % 科特斯系数中的积分表达式 f=1; for i=[0:k-1,k+1:n] f=f.*(t-i); end