贵州师范学院毕业论文（设计）

for j=0:i

x(j+1)=a+j*h; s=1; if j==0 s=s; else for k=1:j s=s*k; end end r=1; if i-j==0 r=r; else

for k=1:(i-j) r=r*k; end end

if mod((i-j),2)==1 q=-(i*s*r); else q=i*s*r; end

17

贵州师范学院毕业论文（设计）

y=1; for k=0:i if k～=j

y=y*(sym('t')-k); end end l=int(y,0,i); C(j+1)=l/q;

z=z+C(j+1)*f1(m,x(j+1)); end

g=(b-a)*z

1）当输入a?0,b?1，n?1,m?2时，即在MATLAB命令窗口输入 >> NCotes(0,1,1,2)即可得用梯形公式的积分值和相应科特斯系数如图1.1

2）当输入a?0,b?1，n?2,m?2时，即在MATLAB命令窗口输入 >> NCotes(0,1,2,2)即可得用辛浦生公式的积分值和相应科特斯系数 如图1.2

3）当输入a?0,b?1，n?4,m?2时，即在MATLAB命令窗口输入 >> NCotes(0,1,4,2)即可得用科特斯公式的积分值和相应科特斯系数 如图1.3

18

贵州师范学院毕业论文（设计）

图 1.1

图 1.2

19

贵州师范学院毕业论文（设计）

图1.3

2 复化梯形求积公式的MATLAB实现

通过f(x)的n?1个等步长节点逼近积分

其中,xk?a?kh,x0?a, 程序二：

function s=trapr1(f,a,b,n) % f是被积函数；

% a,b分别为积分的上下限； % n是子区间的个数； % s是梯形总面积；

20