计算方法-牛顿迭代法的应用

牛顿迭代法的应用一、牛顿法简介牛顿迭代法(Newton'smethod)又称为牛顿-拉夫逊方法(Newton-Raphsonmethod),它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x)=0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根，此时线性收敛，但是可通过一些方法变成超线性收敛。该方法广泛用于计算机编程中。简单迭代法是用直接的方法从原方程中隐含地解出x，从而确定出中3)。而牛顿迭代法是用一种间接而特殊的方法来确定中3)的。下面具体推到牛顿迭代公式。假设气是非线性方程为f(x)=0的一个近似根，把f(x)在气处作泰勒展开:TOC\o"1-5"\h\zf(x)=f(x)+f'(x)(x-x)+f"(x-x)2+kkk2!k若取前两项来近似代替f(x)(称为f(x)的线性化)，则得近似的线性方程f(x)-f(xk)+f)(x-xk)=0设f<xk)丰0，令其解为xk+1，则得x=x-4(1)k+ikf(x)k这称为f(x)=0的牛顿迭代公式。\o"CurrentDocument"f(x)…、一它对应的迭代方程为x=x-£〈显然是f(x)=0的同解方程，故其迭代函数为f(x)中(x)=xk-l^xL(f'(x)古0)k在f(x)=0的根a的某个邻域R(\x—以l<5)内，f(x)R0

If”3)1•If3)1、L<1I中'(x)I=If(x)|2在If”3)1•If3)1、L<1I中'(x)I=If(x)|2二、牛顿法应用1、利用牛顿迭代法计算附息国债的实时收益率附息国债实时收益率的理论计算公式附息国债实时收益率也就是附息国债的内含收益率（或内部收益率，IRR），内含收益率以下面的公式为基础计算回:P=^~^+九.+.』+...+.I+.b(1+y)w(1+yJw+1(1+yJw+2(1+y\+n-1(]+y\+n-1式中：y:附息国债实时收益率（内含收益率）；w：交易结算日至下一利息支付日的天数以上一利息支付日（LIPD）至下一利息支付日（NIPD）之间的天数，即:交易结算日至下一利息支付日的天数w交易结算日至下一利息支付日的天数w=365天（366天）发行一年后的除权价:p=p(1+x)t-11.83发行二年后的除权价:p=p(1+x)-发行一年后的除权价:p=p(1+x)t-11.83发行二年后的除权价:p=p(1+x)-11.83发行九年后的除权价:p=p(1+x)-11.83最后一年的到期本息为：100=p9（1+x）-11.83上面一共10个方程可决定10个未知数：p,p,p,012p9，将以上10个方程整理，消去p1,p2，…,p9得到附息国债实时收益率的实际计算方法下面以696券（696券于1996年6月14日发行，10年期，票面年利率为11.83%，附息国债）为例，介绍附息国债实时收益率的具体计算方法。设1997年6月14日之前的某一天购入696券的价格为P，每百元每年的实时收益率为x，0则:离第一次发放利息（1997年6月14日）的时间为七（年数），

_11.8311.8311.8311.83100则:p0=E+(^*+(^B(^石+(^^•从以上例子可以看出，设距离兑付期为n年(取整数)的附息国债到期本金为M，每年票面利息为：，购入价格为P的国债的到期收益率为x，那么：0y0U+XF+(1+XF•此方程式中已知购入价格P，每年票面利息为c，到期本金M，剩余年数取0整为「(注:这里的「和(5)式中的n有区别)，t为当前交易日至下一附息日的天数除以上一次利息支付日至下次利息支付日的天数，即t为当前日至下一个利息支付日的年数。这是一个非线性代数方程，一般不能直接求解，只能采用逐次通近的迭代方法求解，即当知道解的某一初始近似值X0后，从它开始逐步逼近所要求的解x，我们选用牛顿迭代法解以上方程。利用牛顿迭代法计算牛顿迭代法是迭代法的一种，是求解函数方程的一种有效方法，其基本特征是计算格式简单且收敛较快。给定方程f(x)=0。以及根&的初始近似值七，并假设函数f(x)在％的邻域内导数存在。设x0+Ax云是方程f(x)=0。的根，即f(x0+Ax)=f(&)=0成立，将其在七处按泰勒公式展开得:Ax2f(x)+Axf'(x)+Cf+..•=0.002!f”(x)0略去含Ax2及以后各项得近似等式：f(X0)+Axf'(X0)=0x=x+Ax=x—f(x)0—

f'(x)

0因而得Axx=x+Ax=x—f(x)0—

f'(x)

0重复应用以上做法可得：x.=x—'(X，(n=0,1,2,...).n此即为牛顿迭代公式。

对于方程：p=W一二一+一一，设y=1+x，整理得:0.=0(1+x)+'(1+x)n+tpyn+t-cyn-cyn-icy2一cy-M-c=0.设：f(y)=pyn+t-cyn-cyn-icy2-cy-M-c贝0：f'(y)=(n+1)pyn+t-i-ncyn-i-(n-1)cyn-22cy-c应用牛顿迭代公式：y=y-ZE2，选取误差w,y的初始近似值为y，有:n+1"〃f'(y)cny1y1=y-尝0f'(y0)若|y-y<w，则停止迭代，取y浇y，否则继续迭代110c1y2=y1一f(yy2=y1一f(y1)对于y=y-f(yn1)，如果|y-y>w则继续y=y-f，否则停止nn-1f，(y)'n0-1c»+1»f，(y)n-1n迭代，得到方程的根y牝y，国债收益率x总y-1。n+12、工程应用(1)工程概况在架空输电线路的设计中不同气象条件下架空路线的弧垂、应力和线长占有十分重要的位置，是输电线路力学研究的主要内容。这是因为架空线的弧垂和应力直接影响着线路的正常安全运行，而架空线的微小变化和误差都会引起弧垂和应力相当大的改变。设计弧垂小，架空线的拉应力就大，震动现象加剧，安全系数小，同时杆塔载荷增大因而要求强度提高。设计弧垂过大，满足对地安全距离所需杆塔高度增加，线路投资增大，而且架空线的风摆、舞动和跳跃会造成线路停电事故，若加大塔头尺寸，必然会使投资再度提高。因此，设计合适的弧垂是十分重要的。架空线的线长和弧垂是档距、高差和架空比载、应力的函数。当气象条件发生变化时，这些参数将会发生变化。气温的升降引起架空线的热胀冷缩，是线长、弧垂、应力发生相应变化。大风时的覆冰架空线上比载增加，应力增大由于弹性变形使架空线线长增加。不同气象条件下架空线的各参数之间存在着一定的关系。为了使用方便，常将各种气象条件下的架空线的应力和有关弧垂档距的变化用曲线表示出来，这种曲线称之为应力弧垂曲线，亦称力学特性曲线。架空输电线路应力弧垂曲线为线路杆塔设计、计算、施工提供技术资料，针对传统的绘制方法存在效率低、质量差等缺点，通过对应力和弧垂的计算及过程分析，考虑了控制气象条件的主要因素，根据架空线状态方程的特点，采用牛顿迭代法，应用Matlab语言编制了架空线应力和弧垂的计算程序，依据该程序，输入线路的基本参数和迭代初值即可迅速而准确地计算出应力和弧垂值，结合实例绘制了应力弧垂曲线，结果表明该方法可以极大地提高应力和弧垂的计算精度，有效克服认为误差，简化线路设计工作，提高线路设计效率和质量。2、应力弧垂曲线2.1架空线的状态方程当气温、载荷变化时，支撑与两悬挂点间的架空线应力可按状态方程求解。悬挂点等高时孤立档的状态方程为：63-[§-EE1222—侦E(-t)]62-聆22'2=0(2)212462212241Ey222Ey222令a=-[61-K项(t2-t1)]，g=rr1则式(2)为：623+A622-B=0(3)式中：E为架空线的弹性系数，丫为导线半径，2为档距，t为气温。式(3)为悬挂点等高时孤立档的状态方程。对于计算等高连续档得架空应力时，只要将孤立的档距用代表档距表示，状态方程就具有相同的形式；同样，对于悬挂点不等高时，只要引进一些相应系数，就可以表示与式(3)同样的结构形式。解出状态方程即可求出应力，进而求出弧垂。2.2计算各气象条件下的应力和弧垂求解不含一次项的一元三次方程，本文采用牛顿迭代法。令V=623+A622-B，其导数为：V=3622+2A62，则牛顿迭代式为：给出迭代初值82（。），算出y（0）、y'（0），利用式（4）迭代求出52（i）；反复进行下去，直至I52（〃+1）-52（〃）|<5（5为很小的证书，文中取10-4）为止；最终得出满足方程52即为架空线的应力。由牛顿迭代法求解架空线应力既简单又准确，比采用计算尺或查表法精度高，可是在牛顿迭代法中，有可能由于初值选择不合适而出现迭代不收敛。因此，事先给定初值，进行有限次迭代。如果没有获得结果，就重定初值进行迭代。根据求解出的应力值，利用弧垂计算公式可计算任何档距下标准弧垂的精确值。采用Matlab语言编制求解架空线应力和弧垂的程序，其流程图如图2所示。图2牛顿迭代流程图2.3绘制导线的应力弧垂曲线(1)各气象类型对应的比载及温度如表1所示：Eg.222482口0Eg.222482口0=8-电「12-Ea(t-1)io2482ni10(5)CEg2/2—、其中令A=8i0-驾M-Ea(tn-ti)，10Eg2/2B=^^，则式(5)可简化为:2482n0气彖类犁温度广C)覆冰厚度Jiri)比载<vhihi2)最VI温0■to0gi最低气温U-400gl.4.平凶气温t)0Ri最人顼述10：0曲最大风速30-50io.m1勺过电压15s0与b皿IbJ外过也上10］：50装ID&150£&10?事故0-200卸(2)以控制气象为第一状态，待求状态为第二状态，将第一状态与第二状态所对应的数据分别代入状态方程式：8n03-A8i02-B=0。以平均气温为第一状态，最大覆冰为第二状态为例，运用迭代法得C程序如下：main(){floatl,f,g1=0.03275,g2=0.06544,A,B,t1=-5,t2=-5,c=68.51,x,x1,x0,dx;scanf("%f”,&l);A=c-(73000*g1*g1*l*l)/(24*c*c)-73000*0.0000196*(t2-t1);B=73000*g2*g2*l*1/24;printf("%f,%f”,A,B);x0=1.0;do{dx=(x0*x0*x0-A*x0*x0-B)/(3*x0*x0-2*A*x0);x0=x0-dx;}while(x0!=x1)printf("\nc=%f\n",x0);f=(g2*l*l)/(8*x0);printf("f=%f”,f);}每隔50m,利用以上程序，分别求得表1中九种气象类型下的应力和弧垂，见表2、3、4。表2地线的应力弧垂曲线数据表[神气单供端音％胡以低FEhl均’成EE75【5gi=S3_76Gi=G5.44t=-5I—-Of0f0fflf.-心13-440..11109-610,09270.:il初.(MQ]7IDO岂I.FT1..;*IOSi6L0.S7mI.J].(JSr:<(：.'.C.f：.1S7.812.3QIQS.即§8.81】.鸵说.可1.27I/.U:;X腮：r.tuIUM11(J.8l朋EK:!H8..1I1.w20043.943.T1ia^07L57at572.M5S_512.泊顼48.355.302.&1IO2_53札的B8_513.73成U5il掐T,l.i■I.OIll'71tjiil?5b0.j|n昭32552.QT8.Ha®.LWlog.si7.8®68_51E.213SQR3.3T9.1183.795_UB103-019.H6G_S87.R3IUU.〔母2276.lEiN.ftJIU!IblII.!Mh：i.\l?10.2-n品大IS哧-•1W:J.;.削1.1.1071.04-11.〔订109.61IJ.116L99].：.：,；■■袖U佻都IUIkII'：1<jlIB.li'jliU.[：l1(：.89枷Qi.rlj.a-i11'：'.1家，南胡.弟J'J.7BM0-J.Ji27.1"<;1.史^：：...-!|1UH.L:l■h.522.1.由UF4.-!?11.^4(：l.如U7.骚IU!1.i：l:H.EF.H.2：!JU.71

表3地线的应力弧垂曲线数据表pJftHUli••律心—W.了]t--5jMii：i=S3b94%_ki=J3.17t-lfiffUu.ln=S3B17Z5nf<*ECJt1<*fEO〔T.，r；0.32G1.IS0.17:冷..，］0.2B为.肆514：I(JQ69.020„7305.oy0.65啊.84。.址M就O.M75.43的.3ELMO51.731.707D.^51.11制73,|.：i顷：间29L.iS治,u.-1.网1*1,1*78.4(i17UGft.63u55.53xya77.373.11mu-，:•温颂84..Hli-ft.SSSfi.i'U4.io攵T?:;4'■1Ta'*300K.607o.ii(©a27e.is74.r..|■j.Ul335璐.386,S37Q,206.3861,O7,1773.9®5.9235032・1日土14CiS.44L59g438.4071..II7.11•4U。：；］DO.的05.7(310.32泡m目11.17E7.6Lg.81伽H.02Iii.SBl.<.I'.se.iHII"IM沸u,9：!啊77.9117.&3B3„5416.知58.0717.8503.fra16.4355077.1：^：M.4Jiil.；l以抵对li.=：tl7&HI?f'■?.:;:!GQ0H：!?.，.祁[iU.7£%.1.::」..□36.Ud60.69■JI..J.l匚和荷皿:,u.•伟GO.IS29.7857.Q7郁70Hil.W比载单位：X10-3N/m-mm2；温度单位：。C；应力