数理与信息工程学院
实验报告
课程名称：现代数值分析
实验室：
实验台号：
班级：信息与计算科学121班
姓名：孙佳丽
实验日期：2014年12月03日
实验名称
复化求积公式实验
实验题目
用复化梯形公式、复化辛普生公式、龙贝格公式求下列定积分，要求绝对误差为 ，并将计算结果与精确解进行比较：
（1） ；（2） 。
实验原理
x=linspace(a,b,n+1);%在区间上等分产生n+1个节点
y=feval(f,x);%求f在x节点上的函数值
s=0.5*h*(y(1)+2*sum(y(2:n))+y(n+1));%复化梯形求积公式
复化辛普生msimp.m
functions=msimp(f,a,b,n)
%用复化辛普森公式求积分
n=n+1;
s=msimp(f,2,3,n);
end
n,I,s,error=I-s
龙贝格公式：
formatlong
f=@(x)2*x./(x.^2.-3);
[s,n]=mromb(f,2,3,0.5e-8)
I=log(6)
error=I-s
运行结果：
（1）
复化梯形公式：
复化辛普生公式：
i =
45
I =
龙贝格mromb.m
function[T,n]=mromb(f,a,b,eps)
%用龙贝格公式求积分
%f是被积函数，a,b是积分区间，R返回积分值，eps控制精度，n返回区间等分数
ifnargin<4,eps=1e-6;end
h=b-a;
R(1,1)=(h/2)*(feval(f,a)+feval(f,b));
实验结果分析：
从区间等分上看，复化梯形公式区间等分份数最多，龙贝格公式区间等分份数最少；从误差上看，复化梯形公式误差最大，龙贝格公式误差最小。故复化梯形公式精度最低，龙贝格公式精度最大。
end
err=abs(R(J+1,J+1)-R(J+1,J)) ;
n=2*n;
end
R;
T=R(J+1,J+1);
运行程序：
(1)
复化梯形公式：
formatlong
f=@(x)2/3*x.^3.*exp(x.^2);
n=1;
s=mtrap(f,1,2,n);
I=exp(4);
while(abs(s-I)>0.5e-5)
n=n+1;
s=mtrap(f,1,2,n);
end
n,I,s,error=I-s
复化辛普生公式：
formatlong
f=@(x)2/3*x.^3.*exp(x.^2);
n=1;
s=msimp(f,1,2,n);
I=exp(4);
while(abs(s-I)>0.5e-5)
n=n+1;
s=msimp(f,1,2,n);
%格式：s=msimp(f,a,b,n) f为被积函数，a,b为积分区间的左右端点，n为区间的等分数，s返回积分近似值
h=(b-a)/n;%小区间长度
x=linspace(a,b,2*n+1);%在区间上等分产生n+1个节点
y=feval(f,x);%求f在x节点上的函数值
s=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(n+1));%复化辛普森求积公式
end
n,I,s,error=I-s
龙贝格公式：
formatlong
f=@(x)2/3*x.^3.*exp(x.^2);
[s,n]=mromb(f,1,2,0.5e-5)
I=exp(4)
error=I-s
(2)
复化梯形公式：formFra bibliotektlongf=@(x)2*x./(x.^2.-3);
n=1;
s=mtrap(f,2,3,n);
复化梯形公式
复化辛普生公式
复化梯形公式
实验内容和步骤：
1、编程
2、运行
程序：
复化梯形mtrap.m
functions=mtrap(f,a,b,n)
%用途：用复化梯形公式求积分
%格式：s=mtrap(f,a,b,n) f为被积函数，a,b为积分区间的左右端点，n为区间的等分数，s返回积分近似值
h=(b-a)/n;%区间长度
I=log(6);
while(abs(s-I)>0.5e-8)
n=n+1;
s=mtrap(f,2,3,n);
end
n,I,s,error=I-s
复化辛普生公式：
formatlong
f=@(x)2*x./(x.^2.-3);
n=1;
s=msimp(f,2,3,n);
I=log(6);
while(abs(s-I)>0.5e-8)
54.59815003314423
s =
54.59815464193976
error =
-4.60879553e-06
龙贝格公式：
（2）
复化梯形公式：
复化辛普生公式：
n =
17
I=
1.791759469228055
s =
1.791764166622229
error=
-4.697394174e-09
龙贝格公式：
n=1;J=0;err=1;
while(err>eps)
J=J+1;h=h/2;S=0;
fori=1:n
x=a+h*(2*i-1);
S=S+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*S;
fork=1:J
R(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1);