几种常用数值积分方法的比较

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

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