end
suk; syms t a b fxn2,su=t^2; for u=1:n
su=su*(t-u); end
su; intf=int(su,t,0,n); y=double(intf); RNC= (((b-a)/n)^(n+3))*fxn2*abs(y)/ suk; else
for k=1:n+1 suk=suk*k; end
suk; syms t a b fxn1,su=t; for u=1:n
su=su*(t-u); end
su; intf=int(su,t,0,n); y=double(intf); RNC= (((b-a)/n)^(n+2))*fxn1*abs(y)/ suk;
end
2．分别利用复化梯形公式和复化Simpson公式计算定积分 I??01?exdx
取n?2,4,8,16，精确值为I?4.006994。
复化梯形公式MATLAB程序如下 clear;
Iexact=4.006994; a=0;b=2;
fprintf('\\n n I Error\\n'); n=1; for k=1:4 n=2*n;
h=(b-a)/n;i=1:n+1;
x=a+(i-1)*h;f=sqrt(1+exp(x)); I=travez_v(f,h);
fprintf('%3.0f .5f .5f\\n',n,I,Iexact-I);
7
2end
其中复化梯形求积法的M文件如下 trapez_v.m
function I=trapez_v(f,h)
I=h*[sum(f)-(f(1)+f(length(f)))/2];
复化Simpson求积公式MATLAB程序如下： clear;
Iexact=4.006994; a=0;b=2;
fprintf('\\n n I Error\\n'); n=1;
for k=1:4,n=2*n h=(b-a)/n;i=1:n+1;
x=a+(i-1)*h;f=sqrt(1+exp(x));
I=(h/3)*(f(1)+4*sum(f(2:2:n))+f(n+1)); if n>2
I=I+(h/3)*2*sum(f(3:2:n)); end
fprintf('%3.0f .5f .5f\\n',n,I,Iexact-I); end
3．计算无穷限广义积分I?用-10,10替换积分限，则有
1??edx，其精确值I?1。
???x2??I?1??10?10e?x2dx
利用复化梯形公式，调用trapez_v.m文件，取n=10,20的MATLAB程序如下：
I=1; a=-10;b=10; fprintf('\\n h n I1\\n');
8
n=0;
for k=1:2;n=n+20; h=(b-a)/n;i=1:n+1; x=a+(i-1)*h;f=exp(-x.^2); I=trapez_v(f,h); I1=1/sqrt(pi)*I;
fprintf('%3.1f .0f .8f\\n',h,n,I1); end
4．取精度为??10?8，分别用Tj?1,j?1?Tj?1,j??和Tj?1,j?1?Tj,j??作为计算停止的条件，用Romberg程序计算I??1.501dx，取精度为10-8,并与精确值1?x比较。然后取精度为10-7，用Tj?1,j?1?Tj?1,j??作为计算停止的条件，观察用龙贝格求积公式计算的结果与精确值的绝对误差是否满足10-7。
MATLAB程序如下
F=inline('1./(1+x)'); [RT,R,wugu,h]=romberg(F,0,1.5,1.e-8,13) syms x
fi=int(1/(1+x),x,0,1.5); Fs=double(fi), wR=double(abs(fi-R)), wR1= wR - wugu 其中龙贝格积分的MATLAB程序
function [RT,R,wugu,h]=romberg(fun,a,b, wucha,m) n=1;h=b-a; wugu=1; x=a;k=0; RT=zeros(4,4); RT(1,1)=h*(feval(fun,a)+feval(fun,b))/2; while((wugu>wucha)&(k x=a+h*(2*j-1); s=s+feval(fun,x); end RT(k+1,1)= RT(k,1)/2+h*s; n=2*n; for i=1:k RT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1); end 9 wugu=abs(RT(k+1,k)-RT(k+1,k+1)); end R=RT(k+1,k+1); 四、实验要求： 1．写出梯形公式、Simpson公式、柯特斯公式求数值积分的算法； 2．写出复化梯形公式和复化Simpson公式Romberg方法求数值积分的算法；3．进一步加深对数值积分的理解。 10
