贵州师范学院毕业论文（设计）
h=(b-a)/n; s=0; for k=1:(n-1) x=a+h*k; s=s+feval('f',x); end format long
s=h*(feval('f',a)+feval('f',b))/2+h*s; 先用M文件定义一个名为f.m的函数：
function y=f(x) if x==0 y=1;
else
y=sin(x)/x;
end
若取子区间的个数
在MATLAB命令窗口中输入
>> trapr1('f',0,1,8) 回车得到 如图2.1
21
贵州师范学院毕业论文（设计）
图2.1
3 复化Simpson求积公式的MATLAB实现 程序三：
function s=simpr1(f,a,b,n) % f是被积函数；
% a,b分别为积分的上下限； % n是子区间的个数；
% s是梯形总面积,即所求积分数值； h=(b-a)/(2*n); s1=0; s2=0; for k=1:n
x=a+h*(2*k-1); s1=s1+feval('f',x);
end
22
贵州师范学院毕业论文（设计）
for k=1:(n-1) x=a+h*2*k; s2=s2+feval('f',x);
end
s=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3; 先用M文件定义一个名为f.m的函数：
function y=f(x) if x==0 y=1;
else
y=sin(x)/x;
end 若取子区间个数
时
在MATLAB命令窗口中输入
>> simpr1('f',0,1,8) 回车得到 如图3.1
图3.1
23
贵州师范学院毕业论文（设计）
4 龙贝格积分方法的MATLAB实现 构造数表来逼近积分?af(x)dx?R(J,J)
其中。R(J,J)表示数表的最后一行，最后一列的值。 程序四：
function [R,quad,err,h]=romber(f,a,b,n,delta) % f是被积函数
% a,b分别是积分的上下限 % n+1是T数表的列数 % delta是允许误差 % R是T数表 % quad是所求积分值 M=1; h=b-a; err=1 J=0; R=zeros(4,4);
R(1,1)=h*(feval('f',a)+feval('f',b))/2 while ((err>delta)&(J b24