贵州师范学院毕业论文（设计）
for p=1:M
x=a+h*(2*p-1); s=s+feval('f',x); end
R(J+1,1)=R(J,1)/2+h*s; M=2*M; for K=1:J
R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1); end
err=abs(R(J,J)-R(J+1,K+1));
end
quad=R(J+1,J+1)
先用M文件定义一个名为f.m的函数：
function y=f(x) if x==0 y=1;
else
y=sin(x)/x;
end
在MATLAB命令窗口中输入
>> romber('f',0,1,5,0.5*(10^(-8))) 回车得到 如图4.1
25
贵州师范学院毕业论文（设计）
图4.1
5 高斯-勒让德求积公式的MATLAB实现 程序五：
function [A,x]=Guass1(N) i=N+1;
f=((sym('t'))^2-1)^i; f=diff(f,i); t=solve(f); for j=1:i for k=1:i
X(j,k)=t(k)^(j-1); end
if mod(j,2)==0 B(j)=0;
26
贵州师范学院毕业论文（设计）
else B(j)=2/j; end end X=inv(X); for j=1:i A(j)=0; x(j)=0; for k=1:i
A(j)=A(j)+X(j,k)*B(k); x(j)=x(j)+t(j); end x(j)=x(j)/k; end
function g= GuassLegendre (a,b,n,m)
% a,b分别是积分的上下限； % n+1为节点个数；
% m是调用f1.m中第几个被积函数； [A,x]=Guass1(n); g=0; for i=1:n+1
y(i)=(b-a)/2*x(i)+(a+b)/2; f(i)=f1(m,y(i));
27
贵州师范学院毕业论文（设计）
g=g+(b-a)/2*f(i)*A(i); end
用M文件分别把上面两个自定义函数定义为名为Guass1.m函数和GuassLegendre.m函数
用M文件定义一个名为f1.m的函数function f=f1(i,x) g(1)=sqrt(x); if x==0 g(2)=1; else
g(2)=sin(x)/x; end
g(3)=4/(1+x^2); f=g(i);
在MATLAB命令窗口中输入
>> GuassLegendre (0,1,2,2) >> GuassLegendre (0,1,3,2) 回车得到 如图5.1
28