贵州师范学院毕业论文（设计）
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 贵州师范学院毕业论文（设计） 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 贵州师范学院毕业论文（设计） 图5.1 29