复化辛普森公式Excel-VBA版
复化梯形公式和复化Simpson公式
数值计算方法上机题目3一、计算定积分的近似值：e2= xe x dx要求：（1）若用复化梯形公式和复化Simpson 公式计算，要求误差限= 110-7，分别利用他们的余项估计对每种算法做出步长的事前估计；（2）分别利用复化梯形公式和复化 Simpson 公式计算定积分；（3）将计算结果与精确解比较，并比较两种算法的计算量。
1.复化梯形公式程序：程序 1（求 f（x）的 n 阶导数：syms xf=x*exp（x）%定义函数f（x）n=input（'输入所求导数阶数:'）f2=diff（f,x,n）%求f（x）的n阶导数结果1输入n=2f2 =2*exp（x） + x*exp（x）程序 2：clearsyms x %定义自变量xclcf=inline('x*exp(x)','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('(2*exp(x) + x*exp(x))','x') %定义f(x)的二阶导数,输入程序1里求出的 f2即可。
f3='-(2*exp(x) + x*exp(x))' %因fminbnd()函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=5*10^(-8) %精度要求值a=1 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项if abs(Rn)<e % 用余项进行判断break % 符合要求时结束endendh=(b-a)/n %求hTn1=0for k=1:n-1 %求连加和xk=a+k*h Tn1=Tn1+f(xk)endTn=h/2*((f(a)+2*Tn1+f(b)))z=exp(2)R=Tn-z %求已知值与计算值的差fprintf('用复化梯形算法计算的结果 Tn=')disp(Tn) fprintf('等分数 n=') disp(n) %输出等分数fprintf('已知值与计算值的误差 R=')disp(R)输出结果显示：用复化梯形算法计算的结果Tn= 7.3891等分数n=7019已知值与计算值的误差R= 2.8300e-0082.Simpson 公式程序：程序 1 （求 f（x）的 n 阶导数）：syms xf=x*exp（x）%定义函数f（x）n=input（'输入所求导数阶数:'） f2=diff（f,x,n）%求f（x）的n阶导数结果1输入n=4f2 = 4*exp(x) + x*exp(x)程序 2：clcclearsyms x %定义自变量xf=inline('x*exp(x)','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('(4*exp(x) + x*exp(x))','x') %定义f(x)的四阶导数，输入程序1里求出的f2即可f3='-(4*exp(x) + x*exp(x))'% 因 fminbnd （）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值e=5*10^(-8)% 精度要求值a=1 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) 最大值点对应的x值for n=2:1000000% 求负的四阶导数的最小值点，也就是求四阶导数的% 求等分数 nRn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1) %计算余项if abs(Rn)<e break end end h=(b-a)/n% 用余项进行判断% 符合要求时结束%求 hSn1=0Sn2=0 for k=0:n-1xk=a+k*h%求两组连加和xk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b)) %因Sn2多加了k=0时的值，故减去f(a) z=exp(2)R=Sn-z %求已知值与计算值的差fprintf('用Simpson公式计算的结果 Sn=')disp(Sn)fprintf('等分数 n=')disp(n)fprintf('已知值与计算值的误差 R=')disp(R)输出结果显示：用Simpson 公式计算的结果Sn= 7.3891等分数n=24已知值与计算值的误差R= 2.7284e-008用复化梯形公式计算的结果为：7.3891，与精确解的误差为：2.8300e-008。
复化梯形公式,复化辛普森公式,复化柯特斯公式
复化梯形公式,复化辛普森公式,复化柯特斯公式
复化梯形公式、复化辛普森公式和复化柯特斯公式都是用来计算定积分的近似值的方法。
1. 复化梯形公式：将积分区间分成若干个小区间，在每个小区间上用梯形面积近似代替该小区间的曲边梯形面积，然后将这些梯形面积相加，得到积分的近似值。
2. 复化辛普森公式：将积分区间分成若干个等分小区间，在每个小区间上用矩形面积近似代替该小区间的曲边梯形面积，然后将这些矩形面积相加，得到积分的近似值。
3. 复化柯特斯公式：将积分区间分成若干个等分小区间，在每个小区间上用切线段长度近似代替该小区间的曲边梯形面积，然后将这些切线段长度相加，得到积分的近似值。
这三种方法都是通过将积分区间分成若干个小区间，然后在每个小区间上用近似方法计算该小区间的曲边梯形面积，最后将这些近似值相加得到积分的近似值。
它们的精度和误差都与分区间的大小有关。
高斯-勒让德积分公式在隧道中线坐标计算中的应用
高斯-勒让德积分公式在隧道中线坐标计算中的应用周小利，钟庆丰，刘清政，刘恒杰（中铁工程装备集团技术服务有限公司，河南郑州 450000）［摘要］地铁隧道施工中需要确保建筑限界，尤其是盾构施工测量中需要控制盾构与隧道设计中线的偏差在允许范围内，因此隧道设计中线坐标的计算尤为重要。
文章采用5节点的高斯-勒让德积分公式进行隧道设计中线坐标计算，能达到优于1mm的精度。
通过Excel VBA编程设计了地铁隧道中线坐标计算程序，以某地铁线路的隧道中线坐标计算为例进行验算。
结果表明，该程序计算快速且精确可靠，在工程中有较好的实用性。
［关键词］地铁线路；隧道中线；高斯-勒让德公式；坐标计算［中图分类号］O241.4 ［文献标识码］A ［文章编号］1001-554X（2021）-0083-05 Application of Gauss Legendre quadrature formula in calculation oftunnel centerline coordinatesZHOU Xiao-li，ZHONG Qing-feng，LIU Qing-zheng，LIU Heng-jie地铁隧道施工中建筑结构不能侵入限界，需控制隧道施工中实际线路与设计线路的偏差在允许范围内。
地铁隧道通常采用盾构施工，盾构开挖前需要将隧道设计中线坐标导入测量系统中，指导盾构的掘进方向。
因此，隧道设计中线坐标的计算显得尤为重要，关乎着施工中盾构与设计线路的偏差［1］。
地铁隧道线路设计通常采用直线+缓和曲线+圆曲线+缓和曲线+直线连接方式的线型。
传统的计算线路坐标的方法是根据不同的线型采用不同的坐标计算公式，通过级数展开式来计算缓和曲线的坐标［2］，计算的精度取决于保留的多项式项数，过程较繁琐。
李孟山等采用积分区间n等分，利用复化辛普森公式推导了适用各种线型的计算通用公式［3］，此方法的计算精度取决于n的取值。
李全信采用3～5节点的高斯-勒让德公式进行线路坐标计算，通过数值分析证明了5节点的计算公式适用于各种线型的坐标计算。
公路通用复化辛普森公式匝道点位坐标计算4800源程序
公路通用复化辛普森公式匝道点位坐标计算4800源程序------------------杭浦高速临平互通---------------------本文利用的是计算公路匝道点位坐标的复化辛普森通用公式数学模型,集直线、圆曲线、回旋线通用，占字符内存较小，计算精度不限的程序一、运行变量名称说明：V=1、2分别进入坐标计算、桩号反算K1、K2-------曲线起点、终点里程F0-----------曲线起点方位角R1、R2------曲线起点、终点半径(ρ左-右+，0为直线)X0、Y0-------曲线起点、终点坐标M------------求和累积次数n的2倍(偶数)，精度迭代次数K------------曲线待求点里程BP-----------求点左右偏距(左-右+)ANG-------- -求点的右斜交角X、Y---------曲线求得坐标FW-----------待求点的即时切线方位角XF、YF-------为需求桩号的点坐标DL、K+O、LP分别为桩号误差、求得桩号、左右偏距(左-右+)当曲线的设计半径较小时，为保证点位计算精度，M(即程序中n的2倍)的取值可适当的大些。
M为偶数，直线时M=2即可，经计算M=16即可满足半径为60的小半径曲线精度。
二、曲线计算程序名: Prog "CURVE"Defm 4V"V=1 2"Lbl 0:{KLW}Lbl 4Q"OPT:M0AB1C2D3E4FH5G6I7J8CR9"=0=>Prog "M"△Q=1=>Prog "AB"△Q=2=>Prog "C"△Q=3=>Prog "D"△Q=4=>Prog "E"△Q=5=>Prog "FH"△Q=6=>Prog "G"△Q=7=>Prog "I"△Q=8=>Prog "J"△Q=9=>Prog "CR"△A"K1"B"K2"C"F0"D"R1"E"R2"F"X0"G"Y0"D≠0=>I=1/D:≠=>I=D△E≠0=>J=1/E:≠=>J=E△AbsD+AbsE=0=>M=2:≠=>M=16△V=2=>L=0:W=90△KL"BP"W"ANG"N=0：Z[1]=0：Z[2]=0：Z[3]=0：Z[4]=0 进入坐标迭代计算Lbl2N=N+1：H=2(K-A)/M：R=NH/2+A：R=C+180/π*(I+(J-I)/2(B-A)*(R-A))*(R-A)Int(N/2)=N/2=>Z[1]=Z[1]+cosR：Z[2]=Z[2]+sinR：≠=>Z[3]=Z[3]+cosR：Z[4]=Z[4]+ sinR△N=M=>Goto3：≠=>Goto 2Lbl3X=F+H/6*(cosC+4Z[3]+2Z[1]-cosR)+Lcos(R+W)Y=G+H/6*(sinC+4Z[4]+2Z[2]-sinR)+Lsin(R+W)V=2=>Goto 6△X"X="◢Y"Y="◢R"FW"=R-360Intg(R/360◢Goto 0Lbl 6 进入桩号求算Pol(T"XF"-X,U"YF"-YO=Icos(J-RAbsO≤1e-4=>O"DL"◢K=K+O◢O"LP"=Isin(J-R◢{TU}Goto 6:≠=>K=K+O:L=0:Goto 4三、数据文件：线元要素数据文件每行为一个线元段，逐句执行赋值，直至不满足、运行完成。
复合辛普森公式
实验5 复合辛普森公式李涛 201226100108 计自1201一、实验目的● 用复合辛普森公式计算积分dx x ⎰+482cos 1，使误差不超过-410（注意所给积分特点，做出相应的处理后再计算）二、实验步骤 1.算法原理复合辛普森原理：将区间],[b a 划分为n 等分，在每个子区间[]1,+k k x x 上采用辛普森公式，若记,2121h x x k k +=+则得 ● ∑⎰-===1)()(n k badx x f dx x f I● ).()]()(4)([6121f R x f x f x f h n n k k k k +++=∑-=+记● ∑-=+++=11)]()(4)([6n k k k k n x f x f x f h S● ],)()(2)(4)([6101121∑∑-=-=++++=n k n k k k b f x f x f a f h称为复合辛普森求积公式，其余项为● .),(),()2(180)(11)4(4∑-=+∈-=-=n k k k k k n n x x f h h S I f R ηη于是当],[)(4b a C x f ∈时，与复合梯形公式相似有 ● ),(),()2(180)()4(4b a f h a b S I f R n n ∈--=-=ηη 易知误差阶为4h ，收敛性是显然的，实际上，只要],[)(b a C x f ∈则可得到收敛性，即 ● ⎰=∞→ban n dx x f S )(lim此外，由于n S 中求积公系数均为正数，故知辛普森公式计算稳定。
2.算法步骤复合辛普森：首先将区间],[b a 划分为n 等分，在每个子区间[]1,+k k x x 上采用辛普森公式，若记,2121h x x k k +=+则得∑-=+++=1021)]()(4)([6n k k k k n x f x f x f h S ])()(2)(4)([6101121∑∑-=-=++++=n k n k k k b f x f x f a f h算法过程：这里将辛普森公式写为Sn()函数，然后在Solve()函数里依次计算S1,S2,S4,S6.......当相邻的精度小于eps 时退出循环，则S2n 保存结果。
数值分析复化Simpson积分公式和复化梯形积分公式计算积分的通用程序
数值分析第五次程序作业PB09001057 孙琪【问题】分别编写用复化Simpson 积分公式和复化梯形积分公式计算积分的通用程序；用如上程序计算积分： I (f )=∫sin (x )dx 40取节点x i ， i =0,…,N,N 为2k ,k =0,1,…,12，并分析误差；简单分析你得到的数据。
【复化Simpson 积分公式】Simpson 法则：∫f (x )dx ≈b −a 6[f (a )+4f (a +b 2)+f (b )]b a 使用偶数个子区间上的复合Simpson 法则：设n 是偶数，x i =a +ih , h =b−a n ,(0≤i ≤n) 则有∫f (x )dx =∫f (x )dx +∫f (x )dx +⋯+∫f (x )dx =∑∫f (x )dx x 2i x 2i−2n 2i=1x n x n−2x 4x 2x 2x 0b a 将Simpson 法则应用于每一个区间，得到复合Simpson 法则：∫f (x )dx ≈h 3b a [f (x 0)+2∑f (x 2i−2)n 2i=2+4∑f (x 2i−1)n 2i=1+f (x n )] 公式的误差项为：−1180(b −a )h 4f (4)(δ) 其中δ∈(a,b)【复化梯形积分公式】梯形法则：对两个节点相应的积分法则称为梯形法则：∫f (x )dx ≈b −a 2b a [f (a )+f (b )] 如果划分区间[a,b]为：a =x 0<x 1<⋯<x n =b那么在每个区间上可应用梯形法则，此时节点未必是等距的，由此得到复合梯形法则：∫f (x )dx =∑∫f (x )dx x i x i−1n i=1b a ≈12∑(x i −x i−1)[f (x i−1)+f (x i )]ni=1 对等间距h=(b-a)/n 及节点x i =a +ih ，复合梯形法则具有形式：∫f (x )dx ≈h 2[f (a )+2∑f (a +ih )n−1i=1+f (b )]b a 误差项为：−112(b −a )h 2f ′′(δ)【算法分析】复合Simpson 法则和复合梯形法则的算法上述描述中都已介绍了，在此不多做叙述。
数值分析63 复化求积公式龙贝格求积公式讲解
只增加了一个分点
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改写为
变步长复化辛普森公式计算积分
2. 编写用变步长复化辛普森公式计算积分()ba f x dx ⎰ 的程序。
用上面编写的程序计算下列积分并分析计算结果（1）0cos xdx π⎰ （2）220cos x x dx （3）⎰10dx x程序：function S=bianfuhuasimpson(fx,a,b,eps,M)% 变步长复合simpson 求积公式% 调用方式： S=fuhuasimpson(@fx,a,b,epsilon)% fx -- 求积函数（函数文件）% a, b -- 求积区间% eps -- 计算精度% M--最大允许输出划分数n=1;h=(b-a)/n;T1=h*(feval(fx,a)-feval(fx,b))/2;Hn=h*feval(fx,(a+b)/2);S1=(T1+2*Hn)/3;n=2*n;% 最好与倒数第三行保持一致（变步长）while n<=MT2=(T1+Hn)/2;Hn=0;h=(b-a)/n;for j=1:nx(j)=a+(j-1/2)*h;y(j)=feval(fx,x(j));Hn=Hn+y(j);endHn=h*Hn;S2=(T2+2*Hn)/3;fprintf(' n=%2d S2=%-12.9f S2-S1=%-12.9f\n',n,S2,abs(S2-S1)); if abs(S2-S1)<epsbreak;elseT1=T2;S1=S2;n=2*n;endendS=S2;% 达到下列条件之一，则运算终止:% (1).abs(S2-S1)<eps% (2).下一次的n>M% 输入1：S=bianfuhuasimpson(inline('sqrt(x)*cos(x)'),0,pi,10e-6,2000)% 输入2：S=bianfuhuasimpson(inline('2*x^2*cos(x^2)'),0,sqrt(pi),10e-6,2000) % 输入3：S=bianfuhuasimpson(inline('sqrt(x)'),0,1,10e-6,2000)输出结果：（1）S=bianfuhuasimpson(inline('sqrt(x)*cos(x)'),0,pi,10e-6,2000)n= 2 S2=-0.016369112 S2-S1=0.944423778n= 4 S2=-0.450266122 S2-S1=0.433897010n= 8 S2=-0.669839370 S2-S1=0.219573248n=16 S2=-0.781318443 S2-S1=0.111479074n=32 S2=-0.837710689 S2-S1=0.056392245n=64 S2=-0.866141900 S2-S1=0.028431211n=128 S2=-0.880440980 S2-S1=0.014299080n=256 S2=-0.887620063 S2-S1=0.007179083n=512 S2=-0.891220052 S2-S1=0.003599989n=1024 S2=-0.893023740 S2-S1=0.001803689S =-0.8930（2）S=bianfuhuasimpson(inline('2*x^2*cos(x^2)'),0,sqrt(pi),10e-6,2000)n= 2 S2=1.076354541 S2-S1=2.092222287n= 4 S2=0.039359358 S2-S1=1.036995183n= 8 S2=-0.430456535 S2-S1=0.469815894n=16 S2=-0.662796649 S2-S1=0.232340113n=32 S2=-0.778823323 S2-S1=0.116026674n=64 S2=-0.836827971 S2-S1=0.058004648n=128 S2=-0.865829756 S2-S1=0.029001785n=256 S2=-0.880330615 S2-S1=0.014500859n=512 S2=-0.887581042 S2-S1=0.007250427n=1024 S2=-0.891206256 S2-S1=0.003625214S =-0.8912（3）S=bianfuhuasimpson(inline('sqrt(x)'),0,1,10e-6,2000) n= 2 S2=0.489859598 S2-S1=0.185121744n= 4 S2=0.579745947 S2-S1=0.089886349n= 8 S2=0.623731522 S2-S1=0.043985575n=16 S2=0.645384849 S2-S1=0.021653327 n=32 S2=0.656091436 S2-S1=0.010706587 n=64 S2=0.661402273 S2-S1=0.005310836 n=128 S2=0.664042680 S2-S1=0.002640407 n=256 S2=0.665357576 S2-S1=0.001314896 n=512 S2=0.666013147 S2-S1=0.000655572 n=1024 S2=0.666340270 S2-S1=0.000327122S =0.6663。