龙贝格公式和辛普森公式和复合梯形公式

MATLAB数值分析实验二(复合梯形、辛普森和龙贝格求积,以及二重积分计算等)

佛山科学技术学院实验报告课程名称_______________ 数值分析________________________实验项目_______________ 数值积分____________________专业班级机械工程姓名余红杰学号2111505010 指导教师陈剑成绩日期月日一、实验目的b1、理解如何在计算机上使用数值方法计算定积分 a f ""X的近似值；2、学会复合梯形、复合Simpson和龙贝格求积分公式的编程与应用。

3、探索二重积分.11 f （x, y）dxdy在矩形区域D = ｛（ x, y） | a _ x _ b, c _ y _ d｝的数值D积分方法。

二、实验要求（1）按照题目要求完成实验内容；（2）写出相应的Matlab程序；（3）给出实验结果（可以用表格展示实验结果）；（4）分析和讨论实验结果并提出可能的优化实验。

（5）写出实验报告。

三、实验步骤1、用不同数值方法计算积xln xdx =-- 0 9（1）取不同的步长h，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两公式的精度。

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

2、给出一种求矩形区域上二重积分的复化求积方法，然后计算二重积分..e"y dxdy，其中积分区域D二｛0乞x岂1,0岂y乞1｝。

1.%lnt_t.m复化梯形：function F = Int_t(x1,x2,n)%复化梯形求积公式% x1,x2为积分起点和中点%分为n个区间，没选用步长可以防止区间数为非整数。

%样点矩阵及其函数值：x = lin space(x1,x2 ,n+1);y = f(x);m = len gth(x);%本题中用Matlab计算端点位置函数值为NaN,故化为零: y(1) = 0;y(m) = 0;%算岀区间长度，步长h：h = (x2 -x1)/n;a = [1 2*o nes(1,m-2) 1];%计算估计的积分值：F = h/2*sum(a.*y);%f.mfun cti on y = f(x)y = sqrt(x).*log(x);%run 11.mclc,clear;%分为10个区间，步长0.1的积分值：F = In t_t(0,1,10);F10 = F%分为100个区间F = In t_t(0,1,100);F100 = F%误差计算W10 = abs((-4/9)-F10);W100 = abs((-4/9)-F100);W = [W10 W100]%复化辛普森：%l nt_s.mfun cti on F = In t_s(x1,x2 ,n)%复化梯形求积公式% x1,x2区间，分为n个区间。

复化梯形公式和复化辛普生公式

return result;
}
void simpson::integration()//实现积分
{
cout<<"输入上限和下限";
cin>>b>>a;
cout<<"输入你要使用simposn法则的数目(即等分数)";
cin>>n;
h=(b-a)/n;
sum_even_terms=0.0;
sum_odd_terms=0.0;
for(k=1;k<n;k++)
{
sum_even_terms+=sine(k*h);
}
for(k=0;k<n;k++)
{
sum_odd_terms+=sine((2*k+1)*h/2);
}
integral=(2.0*sum_even_terms+4.0*sum_odd_terms+sine(b)+1)*h/6.0;
《数值分析》实验报告
姓名
学号
日期
2012．11.20
实验室
设备编号
实验题目
用复化梯形公式和复化辛普生公式求I=∫01sinx/xdx
一实验目的
1．了解复化梯形公式和复化辛普生公式。
2.用复化梯形公式和复化辛普生公式求I=∫01sinx/xdx。
二实验内容
算法：复化梯形公式是Tn=∑h/2[f(xi)+ f(xi+1)]=(b-a)/2n[f(a)+2∑f(xi)+f(b)]记子段[xi,xi+1]的中点为xi+1/2,则复化Simpson公式为Sn=∑h/6[f(xi)+4f(xi+1/2)+ f(xi+1)]=b-a/6n[f(a)+4∑f(xi+1/2)+2f(xi)+f(b)]

数值积分：梯形规则

数值积分：梯形规则数值积分：梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式1.问题描述微积分方法求积有很大的局限性，当碰到被积函数很复杂时，找不到相应的原函数。

积分值在几何上可解释为由 x=a,x=b,y=0和y=f(x) 所围成的曲边梯形的面积。

积分计算之所以有困难，就是因为这个曲边梯形有一条边y=f(x)是曲线。

2.理论与方法依据积分中值定理，底为b-a,而高为f(e)的矩形面积恰等于所求曲边梯形的面积I.f(e)称作区间[a,b]上的平均高度。

这样，只要对平均高度f(e)提供一种算法，便相应地获得一种数值求积的算法。

1.梯形规则(Trapezoidal rule)简单选取区间[a ,b]的中点高度作为平均高度。

取h=b-aa0=⌠(a-b)(x-b)/(a-b)dx=(b-a)/2a1=⌠(a-b)(x-a)/(b-a)dx=(b-a)/2得到：2.辛普森规则(Simpson rule)可视作用a , b与c=(a+b)/2三点高度的加权平均值作为平均高度。

3.复合梯形规则(Composite numerical)设将求积区间[a,b]划分为n等份，步长h=(b-a)/2 ,等分点为xi=a+bi , i=0,1,...,n 所谓复化求积法，就是先用低阶求积公式求得每个子段[xi,xi+1]上的积分值,然后再将它们累加求和，用各段积分之和Ii,i=0,1,n-1作为所求积分的近似值。

复化梯形公式：4.复合辛普森规则(Composite Simpson)记子段[xi,xi+1]的中点为则复化公式为复化Simpson公式：5.龙贝格求积公式(Romberg)龙贝格求积公式也称为逐次分半加速法。

它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。

作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度.在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。

龙贝格算法

2
2）理查森外推加速从梯形公式出发，将区间[a,b]逐次二分可以提高求积公式精度，当[a,b] 分为 n 等份时，若记T�� = T ℎ ，当区间[a,b]分为 2n 等份时，则有 T2�� = T
ℎ 2
。再有泰勒公式展开为： T ℎ = �� + ��1 ℎ2 + ��2 ℎ+ ⋯ +
建立一个命名为Romberg.m的function文件：
function[T]=Romberg(f,a,b,e) T=zeros(10,10); T(1,1)=(b-a)/2*(f(a)+f(b)); for k=2:10, sum=0; for i=1:2^(k-2), x=a+(2*i-1)*(b-a)/2^(k-1); sum=sum+f(x); end T(k,1)=T(k-1,1)/2+[(b-a)/2^(k-1)]*sum; for j=2:k, end %第一列用递推梯形公式 %定义龙贝格函数 %定义10阶的零元矩阵
1 3/2 �� 0
��。
算，并取�� , �� ≈ �� ；否则令 k=k+1，转（3）继续计算。 6）下图为我按照自己的算法所设计的示意表：算法设计表：
k 1 2 3 4 …
h b-a (b-a)/2 (b-a)/4 (b-a)/8 …
�� , 1 T(1,1) T(2,1) T(3,1) T(4,1)
�� , 2
�� , 3
�� , 4
�� , 5

数值分析63 复化求积公式龙贝格求积公式讲解

将积分区间 [a, b] 划分为2n等分, 即将每一个区间
[xk, xk+1]经过二等分增加了一个分点
x k ?1/2
?
1 2
(
x
k
?
x k?1)
在每个子区间 [xk, xk+1]上的积分用辛普森公式 , 得
?x k ? 1 xk
f (x)dx ?
h 6
[
f
(
xk
)
?
4
f
(
xk
? 1/2
)
?
f (xk?1)]
?
(b ? a)5 2880n4
f
(4) (? )
?6.3.2 复化辛普森公式
将积分区间 [a, b] 划分为2n等分, 则
? ? ? I ?
b
n?1
f ( x )dx ?
x2k? 2 f ( x )d x
a
k ? 0 x2k
每个子区间 [x2k, x2k+2]上的积分用辛普森公式 , 得
?x2 k ? 2 x2k
称为复化辛普森公式 . 记
? ? h
n?1
n?1
Sn ?
[f 6
(a) ?
4
k?0
f
( x k ? 1/2 ) ?
2
k ?1
f
(xk ) ?
f (b)]
若 f(x)? C 4[a,b], 其求积余项为 h ? b ? a
n
Rn ( f ?
b ? a ( h )4 180 2
f (4) (? ) ?
f (b)]
称为复化梯形公式 . 记
? h
n?1
Tn
?

6b复合求积公式龙贝格算法

步长折半：[xi , xi+1/2] ， [xi +1/2 , xi+1]
n1
xi xi +1/2 xi +1
h T2 n f ( xi ) f ( xi 1 2 ) f ( xi 1 2 ) f ( xi 1 ) i 0 4 n1 h f ( xi ) 2 f ( xi 1 2 ) f ( xi 1 ) i 0 4 h n1 h n1 1 h n1 f ( xi ) f ( xi 1 ) f ( xi 1 2 ) Tn f ( xi 1 2 ) 4 i 0 2 i 0 2 2 i 0 13
1 I T2 n (T2 n Tn ) 3
I Tn 4( I T2n )
3I 4T2n Tn
1 3
验后误差估计式 I T2 n (T2 n Tn )
当
T2n Tn 时，T2n即为所求的近似值。
1 (T2 n Tn ) 3
是T2n 的修正项，它与T2n 之和比T2n 、 Tn更接近与真值，即它是一种补偿。
|| T T2-T|< 2-T1|<
输出T2
16
举例
计算精度满足 | T2n Tn | 107
I [ f ]=0.946083070367
例：用梯形法的递推公式计算定积分解：

1
0
sin( x ) dx , 要求 x
k
0 1 2 3 4 5 6 7 8 9 10
T (k)
梯形法递推公式
1 h n1 1 h n1 T2 n Tn f ( xi 1 2 ) Tn f ( a ih 0.5h) 2 2 i 0 2 2 i 0

龙贝格公式和辛普森公式和复合梯形公式

实验八数值积分信息与计算科学金融崔振威201002034031一、实验目的：1、掌握数据积分算法设计及程序实现二、实验内容：1、p290-1、p301-2三、实验要求：主程序：复合梯形公式：function [I,step,h2] = CombineTraprl(f,a,b,eps)%f 被积函数%a,b 积分上下限%eps 精度%I 积分结果%step 积分的子区间数if(nargin ==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));endendI=I2;step=n;h2=(b-a)/n;function [I,step,h] = IntSimpson(f,a,b,type,eps)%type = 1 辛普森公式%type = 2 辛普森3/8公式%type = 3 复合辛普森公式if(type==3 && nargin==4)eps=1.0e-4; %精度为0.0001endI=0;switch typecase 1,I=((b-a)/6)*(subs(sym(f),findsym(sym(f)),a)+...4*subs(sym(f),findsym(sym(f)),(a+b)/2)+...subs(sym(f),findsym(sym(f)),b));step=1;case 2,I=((b-a)/8)*(subs(sym(f),findsym(sym(f)),a)+...3*subs(sym(f),findsym(sym(f)),(2*a+b)/3)+ ...3*subs(sym(f),findsym(sym(f)),(a+2*b)/3)+subs(sym(f),findsym(sym(f)),b));step=1;case 3,n=2;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/6)*(subs(sym(f),findsym(sym(f)),x)+...4*subs(sym(f),findsym(sym(f)),(x+x1)/2)+...subs(sym(f),findsym(sym(f)),x1));endendI=I2;step=n;endfunction [q,step]=Roberg(f,a,b,eps)%f是被积函数%a积分上限%b积分下限%eps是精度%q输出结果%step循环次数if(nargin==3)eps=1.0e-4;end;M=1;tol=10;k=0;T=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b)); while tol>epsk=k+1;h=h/2;Q=0;for i=1:Mx=a+h*(2*i-1);Q=Q+subs(sym(f),findsym(sym(f)),x);end;T(k+1,1)=T(k,1)/2+h*Q;M=M*2;for j=1:k;T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);end;tol=abs(T(k+1,j+1)-T(k,j));endq=T(k+1,k+1);step=k;p290-11、（i ）用组合梯形公式和M=10求下列每个积分(a)、dx x 1112)1(--⎰+解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('(1+x^2)^(-1)',-1,1)得出结果：q =1.56996299445358s =20h =0.10000000000000 所以dx x 1112)1(--⎰+值约为1.569962994，步长为0.1，M 为20(b)、dx x ⎰+10))2sin(2(解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('(2+sin(2*x^(1/2)))^(-1)',0,1)得出结果：q =0.35117779429220s =19h =0.05263157894737 所以dx x ⎰+10))2sin(2(值约为 0.351177794，步长为0.052631579，M 为19(c)、⎰425.0/x dx解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('1/(x^(1/2))',0.25,4)得出结果：q =3.00216646875717s =46h =0.08152173913043 所以⎰425.0/x dx 值约为 3.002166469，步长为0.081521739，M 为46(d)、dx e x x -⎰402解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('x^2*exp(-1)',0,4)得出结果：q =7.85012162896417s =44h =0.09090909090909所以dx e x x -⎰402值约为7.850121629，步长为0.090909090，M 为44(e)、⎰20)(2dx x xcox解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('2*x*cos(x)',0,2)得出结果：q =0.80323191187607s =36h =0.05555555555556所以⎰20)(2dx x xcox 值约为 0.803231912，步长为0.055555556，M 为36(f)、dx e x x ⎰-π0)2sin(解：在matlab 窗口中输入>> [q,s,h]=CombineTraprl('sin(2*x)*exp(-x)',0,pi)得出结果：q =-6.738211157945265e-018s =1h =3.14159265358979所以dx e x x ⎰-π0)2sin(值约为 -6.7382111578，步长为3.141592653，M 为12、（ii ）用组合辛普森公式和M=10求下列每个积分(a)、dx x 1112)1(--⎰+解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('(1+x^2)^(-1)',-1,1,3)得出结果：I =1.57079538809119step =5s =0.400000000000000所以dx x 1112)1(--⎰+值约为1.570795388，步长为0.4，M 为5(b)、dx x ⎰+10))2sin(2(解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('(2+sin(2*x^(1/2)))^(-1)',0,1,3)得出结果：I =0.35047636367125step =10s =0.10000000000000 所以dx x ⎰+10))2sin(2(值约为0.350476364，步长为0.1，M 为10(c)、⎰425.0/x dx解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('1/(x^(1/2))',0.25,4,3)得出结果：I =3.00030161515673step =14s =0.26785714285714 所以⎰425.0/x dx 值约为3.000301615，步长为0.267857143，M 为14(d)、dx e x x -⎰402解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('x^2*exp(-1)',0,4,3)得出结果：I =7.84809474499077step =4s =1所以dx e x x -⎰402值约为7.848094745，步长为1，M 为4(e)、⎰20)(2dx x xcox解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('2*x*cos(x)',0,2,3)得出结果：I =0.80494830253761step =6s =0.33333333333333所以⎰20)(2dx x xcox 值约为0.804948303，步长为0.333333333，M 为6(f)、dx e x x ⎰-π0)2sin(解：在matlab 窗口中输入>> [I,step,s] = IntSimpson('sin(2*x)*exp(-x)',0,pi,3)得出结果：I =-6.738211157945265e-018step =2s =1.57079632679490所以dx e x x ⎰-π0)2sin(值约为 -6.738211158，步长为1.570796327，M 为2P301-12、（a ）、dx x x ⎰-2024解：在matlab 窗口中输入>> [q,s]=Roberg('sqrt(4*x-x^2)',0,2)得出结果：>> [q,s]=Roberg('sqrt(4*x-x^2)',0,2)q =3.14155917503659s =9所以，该积分值约为3.1415591750（b ）、dx x ⎰+1021/4解：在matlab 窗口中输入>> [q,s]=Roberg('4/(1+x^2)',0,1)得出结果：q =3.14159266527772s =4所以，该积分值约为3.1415926653通过利用龙贝格公式积分求得结果，两个定积分得出结果的积分速度明显不同。

数值分析63 复化求积公式龙贝格求积公式讲解

起来加以考虑 . 注意到每个子区间 [xk, xk+1]经过二分
只增加了一个分点
1 xk?1/ 2 ? 2 ( xk ? xk?1)
设hn=(b? a)/n, xk=a+kh n (k=0,1,? ,n),在[xk, xk+1] 上用梯形公式得
T1 ?
hn 2
?f
(
xk
)
?
f ? ( xk ? 1 )
? ? ? I ?
b
n?1
f ( x )dx ?
xk ?1 f ( x )dx
a
k?0 xk
?h n?1
I
?
6
[
k?0
f
(
xk
)
?
4
f
( x k ? 1/2
)
?
f ( x k?1 )]
? ? h
n?1
n?1
?
[f 6
(a) ?
4
k?0
f
( x k ? 1/2 ) ?
2
k?1
f
(xk ) ?
f (b)]
f ?(??
k)
Rn (
f
)
?
?
b? a 12
h2
f
??(?
),
? ? [a, b], h ? b ? a .
n
可以看出误差是 h2 阶，且由误差公式得到，当
f(x)? C2[a, b] 时，则有
b
? lim
n? ?
Tn
?
f ( x )dx .
a
即复化梯形公式是收敛的. 事实上只要 f(x)? C[a, b], 则可得到收敛，因为只要把 Tn改写为