复平面上超越方程的数值解法及其应用第24卷第2期2008年3月科技通报BULIETIN0FSCIENCEANDTECHNOI0GYV o1.24No.2Mar.2008复平面上超越方程的数值解法及其应用阙仁波,王奎华.(1.浙江大学宁波理工学院土建分院,浙江宁波315100;2.浙江大学岩土工程研究所,杭州310027)摘要:基于幅角原理和闭合曲线积分,结合MA TLAB,提出了复平面上超越方程的两种数值解法,并编制了相应的程序.这两种方法都能将指定区域内的所有解求出,且无需选定初值.数值结果显示了这两种方法的优越性.关键词:幅角原理;闭合曲线积分;MA TLAB;复平面;超越方程中图分类号:O156.6文献标识码:A文章编号:1001—7119(2008)02—0149一O5 NumericalAlgorithmsForTranscendentalEquationsInAComplexPlaneAndItsApplicationsQUERenbo,WANGKuihua(1.DepartmentofCivilEngineeringandArchitecture,Ningboinstituteoftec hnology,ZhejiangUniversity,Ningbo,315100,China;2.InstituteofGeotechnicalEngineering,ZhejiangU niversity,Hangzhou310027,China)Abstract:BasedonArgumentprincipleandcontourintegral,withtheaidofM A TLAB,twonumericalalgorithms combinedwiththecorrespondingproceduresforsolvingtranscendentalequa tionsinacomplexplanearedeveloped. Withtheproposedalgorithms,allsolutionscanbefoundintheappointeddoma iniftheyexist,what‟Smore,thereisnoneedtochooseinitialvalue.Somenumericalexamplesshowtheiradvant ages.Keywords:argumentprinciple;contourintegra1;MA TLAB;complexplane; transcendenta1equations0引言复平面上超越方程的数值求解,是科学研究和工程中经常遇到的问题,同时它也是数值计算中的一个难题,实数域上超越方程的求解方法很多,而复平面上超越方程的求解则相对困难,一般是将其实部和虚部分开,化为二元实数方程组来求解,但有的超越方程很难将其实部和虚部以显式的形式分开,且这种求解方法也不很理想;其它求解方法如牛顿迭代法,Mull—er插值法,蒙特卡洛法和圆盘迭代法等都需要选取初值,且对初值的选取要求很高,而复平面上超越函数实际上相当于二元函数,对其零点的估计是很困难的,给初值的选取带来困难;王则柯采用Kuhn算法L1求解复平面上超越方程L4],取得了较好的效果,但该方法对解的搜索存在随机性且不能保证收敛,因此会出现漏解;收稿日期:2006—1O—O8作者简介:阙仁波(1979一)男,硕士,讲师,主要从事数值计算,振动与波动理论的研究.150第24卷龙云亮等嘲基于幅角原理和闭合曲线积分,将对超越方程的求解转化为对多项式的根的求解,然后利用Kuhn法求解多项式的根,该方法不需要选取初值,且不会漏解,但该方法采用的Ku—hn法,涉及到不动点,单纯同伦算法和拓扑学等方面的知识,概念很难掌握,计算仍有一定复杂性,不便于广大工程技术人员所掌握和应用,因此有必要寻求更为简便的求解方法.许多数学软件如Maple,MA TLAB,Math—ematica,MathCAD等的出现,给使用者带来了极大的方便,许多计算直接调用库函数即可完成,减少了自己编程的繁重任务,其中MA T—LAB以其强大的功能而倍受使用者欢迎.在求解方程方面,这些软件都有着强大的功能,但对超越方程特别是复平面上超越方程的求解问题仍没能很好地解决,对初值的选取要求高,若初值选得不好,就得不到所需要的解, 此外,即使函数有多个解每次也只能给出一个解,但这些函数用于求解多项式的根,都能得到很理想的结果.因此,可利用其优点,避开其不足.鉴于上述情况,本文对龙云亮等_5提出的方法进行改进,将多项式的求解交由MA TLAB中的函数roots或solve去完成,不再采用Kuhn算法;同时提出了另外一种新的求解方法;并对求解区域内存在极点的情况进行了探索.1方法一1.1算法原理设厂(z)为一超越函数,z为复变量,则厂(z)一0即为超越方程,如果f(z)在简单闭曲线C 上解析且不为零,在C的内部也处处解析,,(一1,2,…,押)为厂(z)在C内的所有零点(级的零点算作个零点),则:d一(1)27r厂(z)~令一dz在(2)式中令k一0,由对数留数和幅角原理可知,闭合曲线C围成的区域内,解的个数为:n一dz㈣根据文献(5),求解方程厂(z)=0等价于求解如下的代数方程:z+n1z一+…+a一0(4)其中a(一1,2,…,n)由下面的方程组确定:S1+a1—02+1a1+2a2—0S3+S2a1+Sla2+3a3:0(5)S+S一1a1+Sn--2a2+…+S1a一1+na一0从方程(4)中求解出z(i一1,2,…,n)即为方程厂(z)一0的解.在文献(5)中,采用Kuhn算法求解方程(4),概念和计算复杂,本文对此进行改进,直接调用MA TLAB中的roots或solve函数求解,概念和计算简单,精度高.1.2算法步骤由上述算法原理可知,对于给定的区域,方法一求解超越方程的流程如下图所示:图1流程图一Fig.1flowchartNo.11.3程序实施图2和图3分别表示矩形和圆形区域.对矩形,假定RR,R.R和R.R.,RR分别平行于实轴和虚轴,R,R.处的坐标分别为[R(1),R(2)],JR(3),R(4)];对圆形:圆心0处的坐标为z.,半径为r,圆形区域用[z.r]的表示.为了第2期阙仁波等.复平面上超越方程的数值解法及其应用151 便于积分,先对积分式做如下处理:图2矩形区域图3圆形区域Fig.2rectangulardomainFig.3circulardomainrr对矩形区域:Ig()dz—Ig(z)dz+JcRIR2Ig(z)dz+Ig(z)dz+Ig(z)dz,即将JR2RJR3RJR4R.闭合曲线积分转化为分段积分之和.rr对圆形:Ig(z)dz—Ig(.+)一r2iIg(.+)dO,即将积分转化为周期为J02的函数的积分,积分变量为角度.采用数值积分法如变步长Simpson法求解式(2)和(3)的值.式(5)可写成:SA=B,其中S一100001200021300●●●●●●●●●●●●●●●sn一1sn一2sn一3…n.B一一S1一S2一如●●●一.A—ala2a3●●●n调用MA TLAB中的矩阵左除命令:A:S\B,即可求出A.广Zl1式(4)可写成:E1a…]JZn2I1==0调用roots命令:roots([1a1a2…a])或solve(…“4-a1一4-…4-a一0‟)即可求出解, (=1,2,…,).为了进一步提高精度,可将上面所得到的解作为初值,调用{solve函数对解进一步精确化,只要上面得到的解具备一定的精度,用fsolve函数就能得到更精确的解.2方法二2.1算法原理(1)设_厂()满足方法一中的条件,对图2所示的圆形区域,由(3)式求出,若>0,则方程在该区域内有个解,设为:,,…,,这个解分布的一个特点是:它们与圆心的距离都小于该区域外的解与圆心的距离,即离圆心最近的个解为1,2, (2)(2)MA TLAB中{solve函数求解方程的一般特点是:优先搜索离初值最近的解.(3)结合(1),(2)可知:若以圆心坐标值为初值,采用fsolve函数求解,则可优先求出这个解,避免了初值选取和搜索的随机性.(4)鉴于fsolve函数对每个初值只能求出一个解,而该方法每次都以圆心坐标值为初值,为了不使零点求重(重复求出已求得的解),采用下列方法处理:若已求出r(0<r<)个解1,2,…,,则以fr():一作为下一次求解的函数,显然,只要,.,…,,足够精确,在该区域内函数()就包括且仅包括除1,2,…,,外的解r+1,…,rI‟不会出现求重的现象.(5)若所选区域不是圆形区域,如图3中曲线L所围的区域,但仍满足方法一中的条件,则可通过求解式(3)判断该区域内是否存在图4任意区域解,若存在,做满足方法Fig.4arbitrarydomain 一中的条件且包含该区域的圆形区域,如图3所示,采用上述方法求出圆形区域内的解,再从这些解中找出属于曲线L所围区域内的解,即为所要求的解.计算中,非圆形区域一般选取矩形,三角形区域等,而其相应的圆形区域则选取其外接圆,对矩形区域,若其外接圆中得到的解的实部介于区间[R(1)R(3)]之内且虚部介于区间ER(2)R(4)]之内,则该解即为矩形区域内的解.2.2算法步骤由上述算法原理可知,对于给定的区域,方法二求解超越方程的流程如图5所示:2.3应注意的问题(1)为了提高精度,应对fsolve函数中op—tions的参数进行很好的设置,特别是如下几项:152科技通报第24卷图5流程图二Fig.5flowchartNo.2optimset(…MaxFunEvals‟,MF,‟Maxlter‟,MI,‟TolFun‟,TF,‟TolX‟,TX),MF和MI尽量大,而TF和TX尽量小.(2)为了防止计算误差及该误差在上面函数fr()中传递和放大,一般求解区域不宜过大,大区域可分几次求解.(3)若超越函数表达式中没有虚数单位i的显式时,在选取初值时不要用实数,否则fsolve 不会在复平面上搜索,但同样强行搜索出个解,得出错解而漏解.,(4)为了防止和排除由于误差原因,所得的解不是指定区域内的解,则用fr一()除以(—,)作为下一次求解的函数,同时将总的循环次数增加一次(求出一个解为一次循环),继续搜索,直到求出个属于该区域内的解.3区域内存在极点的情况上面方法一和方法二中都假定区域内没有极点存在,对于存在极点的情况,本文讨论如下简单的形式.如果厂()在简单闭曲线C上解析且不为零,在C内除有限个极点外也处处解析,则:.fd—一㈤其中,P分别为C内零点和极点的总个数,级零点(或极点)算作个零点(或极点).设f(z)可化为:)一,其中),g()之间不可约,F(),g()在闭曲线C上解析且不为零,在C的内部也处处解析,且该函数的极点全由g()确定,而零点则全由F()确定,即g()在F()的零点处解析且不为零,同样F ()在g()的零点处解析且不为零,则:2zriJf铬Fd—(7)()”…2zriJd一(8)g()…一F一㈩厂()().g()…一I』=鬻d(10)故可分别求解F(),g()的零点,所得结果分别为厂()的毒点和极点.4数值试验本文就区域为矩形和圆形两种情况,采用本文的两种算法,编制了通用程序,程序能实现如下计算功能:(1)寻求复平面上超越方程在矩形和圆形区域上的解;(2)求解出指定区域内的解的数目和解(若解存在);(3)指定初始区域和所希望得到的解的最少数目(如在求解偏微分方程时,由于其固有值的多值性,其结果为无穷级数形式,实际计算中不可能取无穷项,但为了达到精度和收敛要求, 要求级数项数不少于N,就得先求出N个固有值),进行求解,对矩形区域还可指定其搜索方向(如沿实轴正向搜索);(4)指定初始区域和所希望得到的解的数目范围(或精确数目),进行求解,对矩形区域还可指定其搜索方向;例1sin()一cos()一0选取矩形区域:[10—0.1200.1搜索得该区域内有4个解,它们是:第2期阙仁波等.复平面上超越方程的数值解法及其应用153 16.4933614313463213.3517687777567619.6349540849362210.21017612418355例2sinz—sin(O.25+0.25i)一0选取初始矩形区域:[一5—151],要求得到4~6个根,且只沿实轴向扩展搜索,扩缩步长为:0.5,搜索得在矩形区域[一6.5—16.51]内有4个解,它们是:一6.03318530717959+0.25000000000000i一3.39159265358979—0.25000000000000i2.89159265358979—0.25000000000000i0.25000000000000+0.2500000000O000i例3一7r++/cos(z)一0选取初始矩形区域:[10—2162],要求至少得到4个根,且实轴向只能沿正向搜索,取实轴向和虚轴向的扩缩步长分别为:2,1,搜索得在区域[10~5225]内有4个解,它们是:20.27307097843463—3.54578255555609i 16.97351747450560+3.36864136733084i 13.94655211299816—3.09348866465839i 10.52169933970673+2.81230084367126i例4一.一0从初始圆形区域[08]出发,要求至少找到6个解,半径扩缩步长为2,搜索得在圆形区域[012]上,该方程有6个解,它们是:7.34327922525260+8.9311638802642li7.34327922525260—8.93116388026421i4.53640365497353—0.00000000000000i1.85718386020784—0.00000000000000i一0.55439223269329+0.61940577182403i一0.55439223269329—0.61940577182403i例5求在圆形区域ll一7r内的解将函数厂()一毫化为/()一的形式,其中F()一sin()一c.s(),g()一1一COS(27r).对F(z),可得在区域ll:丌上有两个解,它们是:0.78539816339745—0.00000000000000i 一2.35619449019218+0.00000000000001i 对g(),可得在区域ll一7c上有14个解, 它们是:0.00000000736029+0.00003907963406i 0.00000000712500—0.00000007000003i 一0.99996639253976—0.00000025439658i1.00001346715059+0.00006066986133i 一1.00000006259024+0.00000000051669i 0.99999995075209—0.00000020561357i 一1.99953736210391+0.00000041868824i 1.99932026164677+0.0006602728290li 一2.00000784021861—0.00000043313944i 2.00001868485329—0.00001437502527i 一2.99903319504505+0.00000137081311i2.99901119152458+0.00028742037317i——3.00220586510496——0.00387496403953i3.00052016151499+0.00037326060752i经分析发现,厂()在该区域内存在2个零点,7个二级极点(=0,±1,±2,±3)计算结果正确.由上面这些例子可看出,该方法能将指定区域内超越方程的解全部求出,且超越函数计算值与精确解(0)的误差都小于10.,充分显示了该方法简便和精确的优越性.参考文献:1-1]王则柯.不动点算法[M].上海:上海科学技术出版社,l987.[2]王则柯,徐森林.逼近零点与计算复杂性理论fJ],中国科学(A),1984(1),8一l5.1-3]王则柯.Kuhn算法的程序实施及数值试验fJ],数值计算与计算机应用,1981(3),175—181.1-4]王则柯.复平面上计算超越函数零点的一种有效算法fJ],中山大学[Ih然科学版],1981(3),15—21.[5]龙云亮,文希理,谢处方.复平面上超越函数零点的数值计算EJ-I,数值计算与计算机应用,2(1994),88 —92.。