高斯数值积分应用

摘要：高斯数值积分方法应用于非圆弧拱坝多拱梁法程序有一定难度，本文针对不同的非圆弧拱特性提出不同积分变量区间的处理方法，将3节点高斯数值积分方法推广应用在5种非圆弧拱(五心拱、抛物线拱、对数螺旋线拱、椭圆拱、双曲线拱)的拱坝多拱梁法程序中，经对比计算，精度较高.

关键词：高斯数值积分方法；积分变量；积分区间；非圆弧拱；多拱梁法

1问题的提出

高斯数值积分方法是一种节点很少、精确度很高的方法，它的特点是节点不等距，计算精度很高，一般利用正交多项式的有关关系式来确定其节点位置和系数.当节点为n时，其代数准确度可达2n-1次.如节点数为3,则求积公式对于任意5(=2×3-1)次多项式都是准确的，这样的精度完全可应用于拱坝程序中的拱段及梁段的计算.笔者曾在圆弧拱多拱梁法程序中，广泛采用3节点高斯数值积分，取得成功[1].但将这一方法推广应用于非圆弧拱多拱梁法程序，却有一定的难度.

由计算数学可知[2]，一般定积分式与高斯积分式的变换形式为：

式中:l=(b-a)/2，为积分限变换系数；n为高斯节点数；ξi为高斯积分节点坐标；gi为与ξi对应的高斯积分系数.

如所周知，拱圈形、载常数计算公式为：

(1)

(2)

式中:S为弧长.

当采用高斯积分公式且为等截面圆拱时，上式相应改为：

(3)

(4)

式中:φ为与弧长S对应的中心角；r为中心半径；E为坝体弹模；Ii,MLi为与高斯节点i对应的截面惯性矩及静定力矩.从上式可见，由于等截面圆拱r为定值，存在简单的dS=rdφ,S=rφ的关系，所以积分变量由S变为φ，积分区间由弧段0～S变为中心角0～φ.在等截面圆拱计算中，许多参数(如坐标及静定力系等)可直接由中心角用显式求得.积分变量由S改为φ后，可大大简化计算.

然而，对于非圆弧拱，问题要复杂得多.因为非圆弧拱的曲率半径和曲率中心处处都在变化，往往不能用简单的显式来表示某一拱段中心角与弧长的关系.并且某些曲线计算弧长也很麻烦，甚至不易用显式求得.如何采用适宜的积分变量和积分区间，是迫切需要解决的问题.

2不同曲线积分变量区间的选取

2.1基本资料5种非圆弧拱示意见图1.几种非圆弧拱曲线方程及有关公式见表1.

图15种非圆弧拱平面示意

表1几种非圆弧拱曲线公式

类型抛物线对数螺旋线椭圆双曲线

曲线方程y=x2/2Rρ=ρ0eaψ

a=cosβ

曲率半径r=(R2+x2)3/2/R2r=Reaψ

r=a2b2x

[(a-y)2/a4+x2/b4]3/2r=a2b2x

[(a+y)2/a4+x2/b4]3/2

弧长S=1/2{xA/cosψA+RLn[(sinψA+1)/cosψA]}S=R/a(eaψA-1)

dSReaψdψ

备注R为拱冠曲率半径β为切线角,R为拱冠曲率半径a为长轴之半,b为短轴之半a为实轴之半,b为虚轴之半

2.2弧的微分公式表1所列弧的微分(dS)算式，除了对数螺旋线稍简单外，其余3种曲线的算式都比较复杂.若从曲线的一般性质来看，当曲线方程y=f(x)时，则弧的微分为：

式中:y′=tgφ，为函数y=f(x)在点x的导数，φ为过点(x,y)的切线与x轴的交角，不难看出，此角与该点的中心角相等.将y′=tgφ代入上式：

所以dS=dx/cosψ(5)

将式(1)作为抛物线、椭圆、双曲线3种曲线弧的微分一般公式，比表1中所列dS算式要简捷得多.式(1)也适用于其他一阶导数存在的任何曲线.

2.3积分变量区间的选取根据各类曲线的性质，选取3种积分变量区间，分述于后.

(1)对于抛物线、椭圆、双曲线3种曲线，采用式(1)所列弧的微分一般公式.此时积分变量为x，积分区间为与S相应的x变化区间.高斯积分时各项均应乘以1/cosφi,φi为与节点i对应的中心角.如A1算式改为:

(6)

(2)对于对数螺旋线，采用表1所列算式：dS=Reaφdφ.此时积分变量为φ，积分区间为与S相应的φ的变化区间，各高斯积分项i均应乘以Reaφi，如A1算式为：

(7)

上式除了高斯积分项乘以eaφi外，其余与圆弧拱相似.

图2五心拱平面示意

(3)对于五心拱，其中弧段为圆弧，算法与等截面圆拱相同，当其边弧段为变截面时，则中心拱弧线为非圆弧曲线，且不便用显式表达.仔细考察该段曲线，发现它与以R3=(RM+RD)/2为半径的圆弧很相近(RM,RD分别为边弧外、内半径)，此圆弧中心O3在边弧起始截面外弧中心O1与内弧中心O2联线中点处，如图2所示.

边弧段计算时，积分变量为φ，积分区间为与边弧S相应的φ的变化区间φ3，如A1算式为：

(8)

上式形式上与圆弧拱一样.

3算例

以上述3种积分变量区间的选取方式，计算各类曲线的半拱弧长，举例于下.

设采用3节点高斯数值积分方法，节点坐标及高斯积分系数列于表2.如以中心角φ为积分变量，以Δφ为积分区间，则与节点i相对应的φi=Δφ(1+ξi)/2.又如以水平坐标x为积分变量，以Δx为积分区间，则与节点i相对应的xi=Δx(1+ξi)/2.

表2节点坐标及高斯积分系数

节点号i节点坐标ξi高斯积分系数gi

1-0.7745966910.555555582

200.888888896

30.745966910.555555582

3.1求抛物线、椭圆、双曲线拱半拱弧长例1：设抛物线拱拱冠曲率半径R=140m，拱端中心角φA=45.32°，拱端坐标xA=141.5726m，求半拱弧长S，可以采用两种方法.一种是按表1所列S的公式直接计算，这是理论积分后的公式，是精确的.另一种是采用高斯数值积分方法，以x为积分变量，xA为积分区间.

(1)按理论公式

S=1/2{xA/cosψA+RLn[(sinψ\-A+1)/coxψA]}

算得S=162.921371m.

(2)按高斯数值积分方法，应有

算得gi/cosφi=2.301610646,S=(141.5726×2.301610646)/2=162.922485m.

该数值积分值与理论计算值162.921371m相比，仅相差0.001114m，相对误差仅为6.8×10-6.

例2：设椭圆拱拱冠曲率半径R=164.51m，拱端中心角φA=51°，坐标xA=141.5843m，长轴之半236.9m，短轴之半197.42m，求半拱弧长.例3：设双曲线拱拱冠曲率半径R=152.71m，拱端中心角φA=41°，坐标xA=141.58421m，实轴之半954.41m，虚轴之半381.77m，求半拱弧长.

上两例因两种曲线无理论积分公式直接用显式计算弧长S，只能与表1中dS算式的高斯数值积分值相比.用dS算式直接数值积分时，积分变量、积分区间与前述方法一样，仍为x及xA，但每一高斯积分项不用除以cosφi，而是乘以dx前的算式等，可见后一算法较繁.两例成果S及比较见表3.

表3椭圆拱、双曲线拱计算成果比较

类别例2椭圆拱例3双曲线拱

原dS算式成果164.143661158.610046

dS=dx/cosψ成果164.143737158.610031

两种算法差值/m0.0000760.000015

相对误差4.63×10-79.457×10-8

由表3可见两种算法成果非常接近.

3.2求对数螺旋线拱半拱弧长例4：设对数螺旋线拱拱冠曲率半径R=154m，拱端中心角φA=46°，切线角β=52°，a=ctgβ=0.781285626,以中心角φ为积分变量，积分区间为拱端中心角φA，求半拱弧长S.

(1)按理论公式计算

S=R/a(eaψA-1)=171.9726625m

(2)按高斯数值积分计算

两种算法成果差值为0.0000061m，相对误差仅为3.547×10-8.

3.3用近似方法求五心拱边例5：设五心拱半拱边弧夹角10°，外半径RM=290m，内半径RD=193.514m，边弧起点拱厚6.466m，拱端厚度8.466m，求半拱边弧长S.

(1)用较精确的计算公式S=φ3×(RA+2×R3)/3.式中3=(RM+RD)/2;φ3为以O3为近似中心的边弧夹角，以弧度计；RA为边弧拱端点与O3联线的长度，见图2.算得RA=241.56728m，φ3=0.206893551弧度，R3=241.757m，从而算出S=50.00488034m.

(2)用近似计算公式S＝φ3×R3=50.01796429m.

两种算法所得的差值为0.013m，相对误差为2.6×10-4.

由以上成果可知，以R3为半径，以φ3为中心角所得的圆弧与边弧的中心弧很近似，因此在高斯数值积分计算时，可以采用较简单的积分变量φ及积分区间φ3.

4结语

综上可知，高斯数值积分应用于非圆弧拱时，需针对非圆弧拱曲线的性质和特点，选取适宜的积分变量与相应的积分区间，这样常可收到事半功倍之效.如对于抛物线拱、椭圆拱、双曲线拱，采用dS=dx/cosφ的通式，既避免了各种曲线的繁复计算，也便于程序规格化.对于对数螺旋线拱，则利用对数函数微分、积分都简单以及极坐标方程弧的微分的特点，积分变量选取φ而不选取x，既大大简化了计算，也很容易利用圆弧拱的算法稍加变换.对于五心拱，在控制误差足够小的前提下，边弧线采用近似圆弧，大大简化了计算.

将3节点高斯数值积分方法推广应用于各种非圆弧拱坝多拱梁法程序，对于提高拱坝程序的计算速度和精度，有很大价值；对于将来推广应用高斯数值积分于各类复杂结构的分析、计算，也有一定的启发和借鉴作用.

参考文献

[1]黎展眉.高斯数值积分方法在多拱梁法程序中的应用[J].砌石坝技术，1986(2).

[2]北京大学，等.计算方法[M].北京：人民教育出版社，1962.