水文计算中伽玛函数的计算方法研究_郭世兴.pdf

伽玛函数作为水文频率及洪水计算中应用广泛的函数,计算方法众多,适用范围和计算精度各不相同。为探索适合水文计算的高精度伽玛函数计算方法,采用斯特林级数及其衍生的近似式等不同伽玛函数渐进展开公式进行计算分析,比较各种方法的适用范围及计算精度。结果表明,分段多项式法截断误差最小,成果精度最高;其次为R a m a n u j a n渐近展开式和s t i r l i n g前四项式。推荐的高精度伽玛函数计算方法可提高水文成果的精度,为各类涉水工程规划、设计中确定工程规模和管理决策提供快速精准的成果。关键词:水文计算;伽玛函数;渐进展开式;斯特林级数中图分类号:TV 1 2 1.1;P 3 3 3 文献标志码:A 文章编号:1 0 0 0-7 7 0 9(2 0 2 3)0 5-0 0 1 5-0 3收稿日期:2 0 2 2-0 5-1 5,修回日期:2 0 2 2-0 8-0 3基金项目:陕西省水利科技项目(2 0 2 1 S L K J-1 7)作者简介:郭世兴(1 9 7 3-),男,正高级工程师,研究方向为水文水资源和水利水电工程规划,E-m a i l:g u o s x 2 0 2 2q q.c o m1 引言概率统计中的皮尔逊型(简称P-型)曲线在工程水文中应用于设计暴雨、洪水、年径流等频率计算,其概率密度函数及分布函数包含伽玛函数,暴雨推求设计洪水的瞬时单位线法中也存在伽玛(简称)函数问题。函数为特殊函数,在数理统计学、物理学、水文水利、通信、金融、农业等众多学科不同实际问题的求解中应用场景广泛1。由于特殊函数的复杂性,应用中常通过级数展开的形式对其近似赋值表示,赋值的误差包含舍入误差和截断误差,前者是由浮点运算的特点造成的,后者则是由所选择的数值近似模型所确定2。特殊函数更精确的赋值和更高效的运算,一直是数学研究的重要方向。特殊函数近似逼近方法有级数展开、渐近展开、连分式展开、P a d 逼近、C h e v b y s h e v逼近、有理函数逼近等方法3,4,其中渐近级数逼近广泛应用于求解微分方程及函数的积分。J a m e s S t i r l i n g3在研究二项分布的正态逼近时提出了伽玛函数的S t i r l i n g f o r m u l a,随 后 衍 生 出S t i r l i n g s e r i e s。L a p l a c e、R a m a n u j a n、B u r n s i d e、G o s p e r、W i n d s c h i t l、N e m-e s、M o r t i c i及陈超平、陈永昌等推导出不种形式的函数展开式,提高其计算精度。文献5 用简单函数对函数1,2区间逼近,误差在0.1 5%。为此,本文通过不同展开式计算值与文献6 进行比较,根据误差范围、收敛情况及效果等方面进行分析,选出最优方法以提高水文成果的计算精度。2 伽玛函数计算方法及成果2.1 伽玛函数自变量取值范围P-型曲线中,频率与(4/C2S)(CS为偏态系数)呈反比。偏态系数C22表示洪水年际变化的不对称度,选取全国7 4个测站和西南“四江”4 8个大河测站分析结果表明7,CS最大值为2.0 5(澜沧江戛旧站),最小值为1.0 4(雅砻江泸宁站);中小测站选取CS变化较大的陕西省分析:渭河支流葫芦河张村驿站CS最大值为4.9 8,汉江安康站CS最小值为1.5 3;小流域中漆水河安头站CS最大值为7.5,岚河六口站CS最小值为1.3 5。故P-型曲线中伽玛函数自变量的取值为0.0 73.7。P-型曲线当CS2时概率密度函数呈“乙”字形,许多干旱、半干旱地区的中小河流洪水,虽然CS2,但经验柱状图仍呈铃形8,所以对自变量小于3.7的取值范围均需精确计算。瞬时单位线(简称I UH)与(n)呈反比(自变量n相当于线性水库个数或调节次数)。当流域汇流时间参数K一定时,n值由小变大,入流过程受调节次数增多,I UH变化由剧烈而趋于平缓,洪峰滞后而减低9。n与流域面积除平均比降(小数计)的幂函数呈反比,通过理想流域1 0(F=3 0 0 k m2、L=3 4.6 4 k m、J=1 5)计算,河流加权平均比降从0.5到1 0 0,n值从3.5 3单调递减到1.9 2。固定河长和比降,流域面积从3 0 0 k m2增 加 到1 0 0 0 k m2,相 同 的 净 雨,n从2.3 8递增到2.7 4。瞬时单位线参数m2=1/n,当m21时取m2=1,故n最小值为1。综上所述,水文分析P-型曲线和I UH中,函数自变量的取值为0.0 73.7。2.2 不同方法的计算结果函数的定义4 为(z)=0e-ttz-1dt(R e(z)0),文献6 给出1,2 区间的七位有效数字的伽玛函数值,其他区间根据伽玛函数递推公式(x+1)=x(x)进行计算,与利用E x c e l中G a mm a L N()函数计算值的相对误差在1 0-8以上。选取具有代表性的几个函数展开式进行分析比较,结果见表1。表1 伽玛函数数值计算方法的结果表T a b.1 R e s u l t s o f a s y m p t o t i c e x p a n s i o n f o r t h e g a mm a f u n c t i o n自变量x图表值s t i r l i n g首四项式s t i r l i n g首两项式s t i r l i n g首两项改进式分段多项式W i n d s c h i t l式R a m a n u j a n式M o r t i c i式C h a o-P i n g C h e n式0.1 5 6.2 2 0 2 7 2 7 3.8 3 2 1 0-1-4.8 0 7 1 0-2-5.3 0 3 1 0-2-3.3 3 4 1 0-88.5 2 9 1 0-2-6.6 4 8 1 0-3-1.2 8 4 1 0-1 5.1 9 1 1 0-10.2 0 4.5 9 0 8 4 3 51.5 3 6 1 0-1-2.6 3 5 1 0-2-3.0 3 6 1 0-2-4.6 1 8 1 0-85.1 2 0 1 0-2-5.4 7 6 1 0-3-1.8 1 4 1 0-11.2 8 3 1 0-10.4 6 1.9 2 5 2 2 6 77.3 0 8 1 0-2-1.5 2 9 1 0-2-1.8 6 5 1 0-22.5 3 1 1 0-83.1 9 6 1 0-2-4.3 9 9 1 0-3-2.3 3 9 1 0-15.8 9 7 1 0-30.5 0 1.7 7 2 4 5 3 87.9 9 0 1 0-3-1.4 5 4 1 0-3-3.4 9 0 1 0-3-4.1 1 3 1 0-86.0 4 6 1 0-3-1.7 6 0 1 0-3-4.5 0 8 1 0-14.1 5 8 1 0-30.7 2 1.2 6 7 4 7 3 15.7 6 2 1 0-3-7.2 4 6 1 0-4-2.6 2 0 1 0-3-2.8 7 2 1 0-84.6 1 1 1 0-3-1.4 9 7 1 0-3-4.9 1 6 1 0-17.7 7 8 1 0-40.9 9 1.0 0 5 8 7 2 01.2 7 4 1 0-38.3 8 7 1 0-4-5.3 5 3 1 0-45.0 3 3 1 0-81.2 6 1 1 0-3-6.6 8 2 1 0-4-7.1 4 9 1 0-11.4 7 7 1 0-41.2 5 0.9 0 6 4 0 2 53.0 3 0 1 0-41.0 2 1 1 0-3-7.1 5 2 1 0-68.9 4 1 1 0-83.5 6 3 1 0-4-2.9 1 3 1 0-4-9.8 7 0 1 0-15.7 3 1 1 0-21.3 6 0.8 9 0 1 8 4 59.7 4 9 1 0-58.9 3 1 1 0-46.5 2 6 1 0-52.5 3 1 1 0-81.3 2 4 1 0-4-4.3 9 9 1 0-3-2.3 3 9 1 0-1-1.5 4 5 1 0-31.4 6 0.8 8 5 6 0 4 36.3 3 3 1 0-58.2 7 2 1 0-46.2 4 3 1 0-5-3.6 6 7 1 0-89.1 4 3 1 0-5-2.6 9 1 1 0-3-3.4 8 1 1 0-15.8 9 7 1 0-31.9 5 0.9 7 9 8 8 0 74.3 6 4 1 0-57.6 9 0 1 0-45.3 7 5 1 0-5-4.1 1 3 1 0-86.6 7 0 1 0-5-1.7 6 0 1 0-3-4.5 0 8 1 0-11.8 5 0 1 0-42.0 0 1.0 0 0 0 0 0 08.5 3 3 1 0-65.3 7 2 1 0-4-5.8 4 6 1 0-68.5 1 0 1 0-81.7 8 7 1 0-5-3.2 6 4 1 0-4-9.4 6 8 1 0-11.3 9 8 1 0-43.2 0 2.4 2 3 9 6 5 48.7 7 8 1 0-92.5 0 7 1 0-4-8.5 7 1 1 0-5-4.6 1 8 1 0-81.6 4 4 1 0-6-2.3 9 7 1 0-5-1.8 1 4 1 0-11.2 5 3 1 0-63.3 8 2.9 1 8 3 1 1 3-9.6 6 8 1 0-82.2 8 7 1 0-4-9.0 2 1 1 0-5-5.5 6 3 1 0-81.2 4 0 1 0-6-1.8 2 9 1 0-56.4 4 2 1 0-47.8 0 3 1 0-74.0 0 6.0 0 0 0 0 0 0-1.4 4 7 1 0-71.7 1 4 1 0-4-9.9 1 3 1 0-50.0 0 0 1 005.7 0 1 1 0-7-8.0 4 9 1 0-62.1 4 5 1 0-41.5 7 6 1 0-7注:相对误差=(图表值计算值)/图表值。(1)s t i r l i n g级数及其首两项近似式。伽玛函数应用最广泛的s t i r l i n g级数形式为:(x+1)2 xxe()x(1+11 2x+12 8 8x2-1 3 95 1 8 4 0 x3-5 7 12 4 8 8 3 2 0 x4+O(x5)(1)式中,x为自变量;e为欧拉数;O为截断误差。S t i r l i n g s e r i e s不收敛,不能像泰勒级数展开一样通过无限制的增加项数来提高精度3,一般取前四项计算,文献1 1 比较后认为增加项数误差反而增大,采用首两项计算。由于(1)=1,根据首 两 项 公 式,本 文 推 导 第 二 项 的 分 母 为1 1.8 4 3x(称为首两项改进式)精度更高。(2)分段多项式1 2(第二公式为欧拉互补定理1 3)。即:(x)1x+c0 x2+c1x3+c2 0 x2 20 x0.5 (2)(x)s i n(x)(1-x)0.5x1 (3)式中,ci(i=0,1,2,2 0)为常数。(3)W i n d s c h i t l渐近展开式1 4。即:(x+1)2 xxe()xxs i n h1x()x/2(4)(4)R a m a n u j a n渐近展开式1 5。即:(x+1)2 xe()xx3+12x2+18x+12 4 0()1/6(5)(5)C.M o r t i c i渐近展开式1 6。即:(x+1)2 ex+1e()x+1/2e x p(11 2x-11 2x2+2 93 6 0 x3-34 0 x4)(6)(6)C h a o-P i n g C h e n渐近展开式1 7。即:(x+1)2 xxe()x1+11 2x3+2 4/7x-1/2()x2+5 3/2 1 0(7)由表1可知,水文计算中函数的图像在第一象限,且在(0,1.4 6 区间单调递减,1.4 6,)区间单调递增;分析图表值与不同函数式计算61水 电 能 源 科 学 2 0 2 3年 第4 1卷第5期郭世兴:水文计算中伽玛函数的计算方法研究值的相对误差值,得出分段多项式计算精度最好(相对误差在1 0-8),其次n1采用R a m a n u j a n式(相对误差在1 0-3),n1采用s t i r l i n g前四项式(相对误差在1 0-4),在n1.3 6时s t i r l i n g首两项改进式略优于s t i r l i n g前四项式。当n趋近无穷大时,各种方法计算的截断误差均在减小,计算值趋于相同。3 结论a.P-型密度曲线形状的参数CS在洪峰流量计算时取值1.0 4,7.5,对应伽玛函数自变量的取值0.0 7,3.7。计算的()值偏小,影响成果超越 频 率 值p偏 大;I UH适 用 于 流 域 面 积3 0 01 0 0 0 k m2无资料地区的设计暴雨推求设计洪水计算,线性水库个数n与流域面积除平均比降(小数计)的幂函数呈反比,取值不小于1.0。计算的(n)偏小,影响计算的洪峰流量偏大。b.水文计算中函数自变量均大于0,图像在第一象限0,1.4 6 区间单调递减,1.4 6,)区间单调递增。采用分段多项式计算截断误差最小,精 度 最 高;其 次 为0 n 1.0采 用R a m a n u j a n渐近展开式,n1.0采用s t i r l i n g前四项式法(n1.3 6时s t i r l i n g首两项改进式精度高于前四项式法)。参考文献:1 刘颖.不完全伽玛函数计算算法改进与应用D.成都:电子科技大学,2 0 1 5.2 常晓阳.几类特殊函数的赋值分析研究D.上海:华东师范大学,2 0 1 8.3 王竹溪,郭敦仁.特殊函数概论M.北京:北京大学出版社,2 0 0 0.4 HA R R Y B A T EMAN.H i g h e r T r a n s c e n d e n t a l F u n c t i o n s,V o l u m e I I IM.N e w Y o r k:M c G r a w-H i l l B o o k C o m p a n y,I n c.,1 9 5 3.5 赵晓慎,马建琴.皮尔逊型分布逼近、递推和迭代计算J.水电能源科学,2 0 0 6,2 4(4):4 4-4 7.6 谭维炎,张维然.水文统计常用图表M.北京:水利出版社,1 9 8 2.7 孙济良.论水文频率线型选优及参数估计J.水力发电,1 9 9 2(3):1 4-2 0.8 中华人民共和国水利部.水利水电工程设计洪水计算规范:S L 4 4-2 0 0 6S.北京:中国水利水电出版社,2 0 0 6.9 长江水利委员会.水文预报方法:第二版M.北京:水利水电出版社,1 9 9 3.1 0长江水利委员会水文局.水利水电工程设计洪水计算手册M.北京:中国水利水电出版社,1 9 9 5.1 1陈永昌.瞬时单位线法中时段单位线的计算式J.水文,1 9 9 4,1 4(1):4 4-4 7.1 2王光生,宁方贵,肖飞,等.实用水文预报方法M.北京:中国水利水电出版社,2 0 0 8.1 3埃伯哈德蔡德勒,等编.数学指南实用数学手册M.李文林,等译.北京:科学出版社,2 0 1 2.1 4M A B R AMOW I T Z,I A S T E GUN.H a n d b o o k o f M a t h m a t i c a l:F u n c t i o n s w i t h F o r m u l a s,G r a p h s,a n d M a t h e m a t i c a l T a b l e sM.N e w Y o r k:D o v e r P u b l i-c a t i o n s,1 9 7 2.1 5R AMANU J AN S.T h e L o s t N o t e b o o k a n d O t h e r U n p u b l i s h e d P a p e r sM.N a r o s a P u b l i s h i n g H o u s e,N e w D e l h i,1 9 8 8.1 6C R I S T I N E L MOR T I C I.N e w i m p r o v e m e n t s o f t h e s t i r l i n g f o r m u l aJ.A p p l i e d m a t h e m a t i c s a n d c o m-p u t a t i o n,2 0 1 0,2 1 7:6 9 9-7 0 4.1 7 CHAO-P I N G CHE N.A m o r e a c c u r a t e a p p r o x i m a-t i o n f o r t h e g a mm a f u n c t i o nJ.J o u r n a l o f n u m b e r t h e o r y,2 0 1 6,1 6 4:4 1 7-4 2 8.S t u d y o n C a l c u l a t i o n M e t h o d f o r G a m m a F u n c t i o n i n H y d r o l o g i c a l C o m p u t a t i o nGUO S h i-x i n g(S h a a n x i P r o v i n c a l I n s t i t u t e o f W a t e r R e s o u r c e s a n d E l e c t r i c P o w e r I n v e s t i g a t i o n a n d D e s i g n,X ia n 7 1 0 0 0 1,C h i n a)A b s t r a c t:A s a w i d e l y u s e d f u n c t i o n i n h y d r o l o g i c f r e q u e n c y a n d f l o o d c o m p u t a t i o n,t h e G a mm a f u n c t i o n h a s m a n y c a l c u l a t i o n m e t h o d s,a n d t h e s c o p e o f a p p l i c a t i o n a n d c o m p u t a t i o n a c c u r a c y a l s o v a r y.I n o r d e r t o e x p l o r e a h i g h-p r e c i s i o n G a mm a f u n c t i o n c o m p u t a t i o n m e t h o d s u i t a b l e f o r h y d r o l o g i c a l c o m p u t a t i o n,d i f f e r e n t G a mm a f u n c t i o n a s y m p t o t i c e x p a n-s i o n f o r m u l a s s u c h a s S t i r l i n g s e r i e s a n d i t s d e r i v e d a p p r o x i m a t i o n f o r m u l a a r e u s e d f o r c o m p u t a t i o n a n d a n a l y s i s.C o m p a-r i n g t h e a c c u r a c y a n d a p p l i c a t i o n r a n g e o f v a r i o u s m e t h o d s,t h e r e s u l t s s h o w t h a t t h e p i e c e w i s e p o l y n o m i a l m e t h o d h a s t h e s m a l l e s t t r u n c a t i o n e r r o r a n d t h e h i g h e s t a c c u r a c y;T h e s e c o n d i s t h e R a m a n u j a n a s y m p t o t i c e x p a n s i o n a n d t h e S t i r l i n g f i r s t t e t r a n o m i a l.T h e r e c o mm e n d a t o r y h i g h-p r e c i s i o n G a mm a f u n c t i o n c o m p u t a t i o n m e t h o d c a n e f f e c t i v e l y i m p r o v e t h e c o m p u t a t i o n a c c u r a c y o f h y d r o l o g i c a l r e s u l t s,w h i c h c a n p r o v i d e f a s t a n d a c c u r a t e r e s u l t s f o r t h e p l a n n i n g a n d d e s i g n o f v a r i o u s h y d r o p r o j e c t s t o d e t e r m i n e p r o j e c t s c a l e a n d m a n a g e m e n t d e c i s i o n s.K e y w o r d s:h y d r o l o g i c a l c o m p u t a t i o n;G a mm a f u n c t i o n;a s y m p t o t i c e x p a n s i o n;S t i r l i n g s e r i e s71