数值分析实验报告2

实验报告
一、实验名称
复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式及自适应辛普森积分。

二、实验目的及要求
1. 掌握复合梯形求积计算积分、复合辛普森求积计算积分、龙贝格求积计算积分和自适应辛普森积分的基本思路和步骤.
2. 培养Matlab 编程与上机调试能力. 三、实验环境
计算机，MATLAB 软件 四、实验内容
1.用不同数值方法计算积分9
4
ln 1
0-=⎰
xdx x 。

（1）取不同的步长h 。

分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确指比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善。

（2）用龙贝格求积计算完成问题（1）。

（3）用自适应辛普森积分，使其精度达到10-4。

五、算法描述及实验步骤
1.复合梯形公式
将区间[a,b]划分为n 等份，分点x k =a+ah,h=(b-a)/h,k=0,1,...,n ，在每个子区间[x k ,x k +1](k=0,1,...,n-1)上采用梯形公式（），得
)]()([2
)(b f a f a b dx x f b a
+-≈⎰ （） )]()(2)([2)]()([21
1
110b f x f b f h
x f x f h T n k k k n k k n ++=+=∑∑-=+-= （）
),(),(12
)('
'2b a f h a b f R n ∈--
=ηη （） 其中Tn 称为复合梯形公式，Rn 为复合梯形公式的余项。

2.复合辛普森求积公式
将区间[a,b]划分为n 等份，在每个子区间[x k ,x k +1](k=0,1,...,n-1)上采用辛普森公式（），得
)]()2
(4)([6b f b a f a f a b S +++-= （）
)]()(2)(4)([6)]
()()([61
1
102/112/11
b f x f x f b f h
x f x f x f h S n k k n k k k k n k k n +++=++=∑∑∑-=-=+++-= （） ),(),()2
(180)()
4(4b a f h a b f R n ∈-=
ηη （） 其中Sn 称为复合辛普森求积公式，Rn 为复合辛普森求积公式的余项。

3.龙贝格算法 统一的公式：
)(1
41
)2(144)(11h T h T h T m m m m m m -----= （）
经过m （m=1,2...）次加速后，余项便取下列形式：
...)()2(22)1(21+++=++m m m h h I h T ξξ （）
上述处理方法通常称为理查森外推加速法。

设以)(0k T 表示二分k 次后求得的梯形值，且以)
(k m T 表示序列｛)(0k T ｝的m 次加
速值，则依递推公式（）可得
,...2,1,1
41144)()
(1)1(1)(=---=-+-k T T h T
k m m
k m m m k m
（） 公式（）也称为龙贝格求积算法，计算过程如下：
（1）取k=0,h=b-a,求)]
()([2)0(0b f a f h
T +=。

令k →1（k 记区间[a,b]的二分次数）。

（2）求梯形值T 0((b-a)/2k ),即按递推公式（）计算)(0k T 。

∑-=++=1
2/12)(221n k k n n x f h T T （）
（3）求加速值，按公式（）逐个求出T 值。

（4）若ε<--)
0(1)0(k k T T （预先给定的精度），则终止计算，并取I T k ≈)0(；否则
令k k →+1转（2）继续计算。

4.自适应积分方法
设给定精度要求0>ε，计算积分dx x f f I b
a
⎰=)()(的近似值。

先取步长h=b-a ，
应用辛普森公式有
),(),()2
(180),()()(4
4b a f h a b b a S dx x f f I b
a ∈--
==⎰ηη （）
表区间[a,b]对分，步长h 2=h/2=(b-a)/2，在每个小区间上用辛普森公式，得
),(),()2
(180),()(4
422b a f h a b b a S f I ∈--
=ξξ （） 上式即为
),(),()4
(180),()(4
42b a f h a b b a S f I ∈--
=ξξ （） 将（）与（）比较得 212215
1),(),(151),()(S S b a S b a S b a S f I -=-≈
- 则期望得到
ε<-),()(2b a S f I （） 此时可取S 2(a,b)作为dx x f f I b
a ⎰=)()(的近视，则可达到给定的误差精度ε。

如果不行，则细分区间，进行计算。

六、调试过程及实验结果
取不同的步长，得到的不同结果如下表：
七、总结
通过本次学习Matlab ，掌握了复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式及自适应辛普森积分的程序和算法，为以后处理数据提供一种更加简便，准确的方法。

八、附录（源程序清单）
1.复合梯形
function s=fuhetixing(f,a,b,n) %f 为被积分函数 %a ，b 是积分上下限 %n 是子区间个数 %s 是积分值 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;
2.复合辛普森
function S=Comsimpson(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
for k=1:(n-1)
x=a+h*2*k;
s2=s2+feval('f',x);
end
format long
S=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3;
3.龙贝格
function [T,quad,err,h]=Romberg(f,a,b,n,delta) %f为被积分函数
%a，b是积分上下限
%n+1是T数表的列数
%T表示T数表
%quad是所求积分值
%delta是设定的允许误差限
m=1;
h=b-a;
err=1;J=0;
T=zeros(n,n);%定义T表初始值
T(1,1)=h*(feval('f',a)+feval('f',b))/2;
while ((err>delta)&(J<n))
J=J+1;
h=h/2;
s=0;
for k=1:m
x=a+h*(2*k-1);
s=s+feval('f',x);
end
T(J+1,1)=T(J,1)/2+h*s;
m=2*m;
for i=1:J
T(J+1,i+1)=T(J+1,i)+(T(J+1,i)-T(J,i))/(4^i-1);
end
err=abs(T(J,J)-T(J+1,i+1));
end
format long
quad=T(J+1,J+1)
err
T
4.自适应辛普森求积公式
function s=S_Adapt_Simpson(a,b,err,M)
%input: a--下限
% b--下限
% err--the tolerance（容差）
% m--初始设置的步数
format long
h=(b-a)/M;%步距
s=0;
for i=1:M
x=a+(i-1)*h;
y=a+i*h;
e=abs(simpson(x,y,2)+simpson(x,y,1))/10;
j=1;
while(e>=err) %循环直到to<tol为止
j=j+1;
e=(abs(simpson(x,y,2^j)-simpson(x,y,1)))/10; %精度测试式
end
s=s+simpson(x,y,2^j);
end。