化工数值计算-2

1、非线性代数方程的数值解用代数方法求解一元线性方程式或一元二次方程式是比较简单的,但求解一元三次或高次的代数方程式却相当麻烦。这类问题在化学化工计算中并不少见。例如,适用于中压范围的范德华方程式为:如果要计算在温度T和压力P时,某种气体的摩尔体积V,将会遇到如下形式的一元三次方程式:再如计算一定温度和压力下,气体反应达平衡时各组分的含量,也时常要解一元三次、一元四次甚至更高次的代数方程。这些高次方程,求解析解不可能,只能用数值方法求近似解。以计算机为工具,利用数值计算法处理以上问题,很容易得到极为准确的近似解。这些非线性代数方程的数值计算方法及计算程序,即为本章的主要内容。2.1逐步扫描法求根的

2、初始近似值通常把方程f(x)=0的解叫做方程的根,也称为函数f(x)的零点。若f(x)是n次多项式,对应的方程为n次代数方程,这时方程的根也称为多项式的根。根有实根和复根之分,这里仅研究实根的求法。可分两步来求方程的根,先找出根的某个粗略近似值,又称之为“初始近似值”,进一步将初始近似值逐步加工成满足精度要求的结果。下面先讨论求根的初始近似值的方法。2.1.1方法概述设待解方程为:在直角坐标系中给出相应于y=f(x)的曲线,如图2-1所示。显然,此曲线与x轴的交点就是方程f(x)=0的根。若函数f(x)在区间(a,b)连续,且f(a)与f(b)异号,则区间(a,b)内必定至少有一个实根。若函数

3、f(x)在区间(a,b)连续并单调(单调上升或单调下降),则在区间(a,b)必定只有一个实根(见图2-2)。这时,选定一个步长h,并计算函数f(a),f(a+h),f(a+2h),的值,直至相邻两个函数值异号,则所求的根必在这两个x值之间。可取xi或xi+h作为所求根的近似值。这个方法叫迈步法或逐步扫描法。图2-1y=f(x)曲线图2-2单调变化的y=f(x)曲线2.1.2程序框图上述算法的程序框图如图2-3所示。图2-3逐步扫描法求根近似值的程序框图程序框图中的主要变量:A初值H步长X1区间始点x1X2区间终点x2Y1区间始点函数值f(x1)Y2区间终点函数值f(x2)计算过程从初值A开始,

4、每增加一个步长H,便对该区间的始点和终点的函数值符号作一次检验。若两函数值乘积大于零,说明该区间内无实根。这时,把区间的终点作为下一次迈步的始点。如此反复向前迈步,直至两函数值的乘积小于或等于零,并把两函数乘积小于或等于零的区间的始点作为方程式根的近似值。显然,步长缩小,精度提高。因此,当精度要求较高时,步长要求很小,计算机循环计算的次数将会很大,因此,用这种方法求根是不切实际,它只能用于求根的初始近似值。2.1.3计算实例例2-1用迈步法求方程x3-x-1=0的根的近似值。解选初值A=0,步长H=0.1源程序:*Example2-1-Eg2-1.frm*DefDblA-H,O-ZPrivat

5、eSubCommand1_Click()ClsA=0:H=0.1PrintPrintA=;A:PrintH=;HX1=A:Y1=Y(X1)OKstop=1DoWhileOKstop0X2=X1+H:Y2=Y(X2)OKstop=Y1*Y2IfOKstop0ThenX1=X2:Y1=Y2EndIfLoopPrintX=;X1EndSub*FunctionY(X)*Y=X3-X-1EndFunction执行结果:A=0H=0.1X=1.30源程序中将待求解方程写成y=f(x)=0的形式,函数y=f(x)的计算采用函数子程序,这样当求解方程改变时,只要改变函数子程序即可。2.2求根的精确解用迈步法求

6、得非线性代数方程根的初始近似值后,可用以下介绍的二分法、迭代法、牛顿法、弦截法等寻求根的精确值。2.2.1二分法2.2.1.1方法概述设求解方程f(x)=0通过逐步扫描法已知有根区间为(x1,x2),如图2-4所示。图2-4二分法求根的精确值现取x1与x2的中点x0,即从而将区间(x1,x2)分为相等的两部分。然后检查f(x0)与f(x1)的符号是否相同,如为同号,根必在x0与x2之间。令如为异号,则根必在x1与x0之间。令图2-5二分法通用计算程序框图然后再取新区间(x1,x2)的中点,重复以上步骤,直至x1与x2之间的距离小于某指定值E为止。取前一区间的中点x0作为方程式根的精确值。这个方

7、法称为“二分法”或“对分法”。若将二分过程无限地进行下去,上述有根区间最终必会收缩于一点x*,x*就是所求的根。用这种方法,可以获得一个近似根的序列:该近似根序列是以x*为极限的。实际计算不可能是无限过程,由于所以只要有根区间(x1,k+1,x2,k+1)的长度小于指定的误差值E,那么就可以认为结果x0,k是能满足方程式f(x)=0的。只要区间(x1,x2)内有根,用二分法一定能够求出结果,因而它是最保险的一种方法,但其最大的缺点是收敛速度慢。2.2.1.2程序框图图2-5是二分法求解非线性代数方程式的通用计算程序框图。这个计算程序包括逐步扫描法求根的初始近似值和二分法求根的精确值。程序框图中

8、的主要变量:E允许误差X0区间中点x0Y0区间中点函数值f(x0)其余符号与迈步法相同。2.2.1.3计算实例例2-2用二分法求方程x3-x-1=0的根。允许误差E=0.0001。解按程序框图(见图2-5)编写源程序。同样,源程序中把函数值的计算作为一个函数子程序来调用。源程序:*Example2-2-Eg2-2.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()ClsA=0:H=0.1:E=0.0001Print:PrintA=;A:PrintH=;H:PrintE=;EX1=A:Y1=Y(X1)-OKstop=1DoWhileOKstop0X2=X1+H

9、:Y2=Y(X2)OKstop=Y1*Y2IfOKstop0ThenX1=X2:Y1=Y2EndIfLoop-OKstop=1DoWhileOKstopEX0=(X1+X2)/2:Y0=Y(X0)IfY1*Y00ThenX1=X0ElseX2=X0EndIfOKstop=Abs(X2-X1)Loop-PrintX=;:PrintFormat$(X1,#.#)EndSub*FunctionY(X)*Y=X3-X-1EndFunction执行结果:A=0H=0.1E=0.0001X=1.3247用上述方法求根,只能得到方程式的一个实根。实际上,若方程式有数个实根时,只需增加一个终值B,并对图2-5

10、作如下修改:当求出第一个实根之后,检查有根区间的终点x2是否小于终值B。若小于终值B,则再施行迈步法求下一个实根的近似值,进而求出精确值。如此反复进行下去,直至达到或超过终值B为止。图2-6是这种算法的程序框图,并用例2-3进行说明。例2-3在理想搅拌釜多釜串联反应器中进行如下液相反应:AR。已知该反应的动力学方程为:rA=kc1.5。其中rA为反应物的消耗速率kmol/(m2min),k为反应速率常数(m3)0.5/(kmol0.5min),c为反应物A的瞬时浓度kmol/m3。若单釜有效容积VR=1m3,各釜容积和温度相同,在该温度下反应速率常数k=0.00138(m3)0.5/(kmol

11、0.5min),物料的体积流量F=3.510-3m3/min,反应物A的起始浓度c0=3kmol/m3,希望A的转化率xp不低于85%,试计算所需的串联釜数和每釜出口反应物浓度和转化率。解根据对理想搅拌釜多釜串联反应器作物料衡算可得:将=VR/F,ri=k代入式(2-3)得:第i釜的转化率为用二分法求式(2-4)的解,然后用式(2-5)计算转化率。在求方程式(2-4)的根时,用前一釜的出口浓度ci-1作为初值。因此,用逐步扫描法求根的近似值时,是沿着反应物浓度由大至小的方向进行迈步的。求出方程式(2-4)的根之后,立即计算转化率,若转化率低于85%,即把求出之根作为下一次求根的初值,再重复上述

12、步骤,直至转化率等于或大于85%为止。其计算程序框图如图2-7所示。程序框图中的主要变量:K反应速率常数kV单釜有效容积VRm3C0反应物A的起始浓度c0kmol/m3CC第i釜反应物A的进口浓度ci-1XP转化率xpI串联釜数F物料体积流量m3/minXI转化率xiT单釜平均停留时间X1、X2、X0、Y1、Y2、Y0、E、H等与图2-5相同。源程序:*Example2-3-Eg2-3.frm*DefDblA-H,O-ZDimVK,T,CCPrivateSubCommand1_Click()ClsVK=0.00138:C0=3:XP=0.85:F=0.0035V=1H=0.2:E=0.001P

13、rint:PrintTab(5);*DATA*PrintK=;VK:PrintC0=;C0PrintXP=;XPPrintF=;F:PrintV=;VPrintH=;HPrintE=;EPrintTab(4);*Result*PrintTab(0);No.;Tab(10);C;Tab(24);XT=V/F:CC=C0:I=1XI=0DoWhileXI0OKstop=Y1*Y2X1=X2-H:C=X1:Y1=Y(C)OKstop=Y1*Y2IfOKstop0ThenX2=X1:Y2=Y1EndIfLoop-OKstop=1DoWhileOKstopEX0=(X1+X2)/2:C=X0:Y0=Y(

14、C)IfY1*Y00ThenX1=X0ElseX2=X0EndIfOKstop=Abs(X2-X1)Loop-XI=(C0-X0)/C0PrintFormat$(I,#);PrintFormat$(X0,0.0000);PrintFormat$(XI,0.0000)IfXI0X2=X1+H:T=X2:Y2=Y(T)OKstop=Y1*Y2IfOKstop0ThenX1=X2:Y1=Y2EndIfLoop-PrintNo.T/KOKstop=100*E:I=0DoWhileOKstopEX0=(X1+X2)/2:T=X0:Y0=Y(T)IfY1*Y00ThenX1=X0ElseX2=X0EndI

15、fOKstop=Abs(X2-X1):I=I+1PrintFormat$(I,#.#);PrintFormat$(X0,#.#)Loop-EndSub*FunctionY(T)*YM=0ForI=1ToNYM=YM+X(I)*Exp(A(I)-(B(I)/(T+C(I)NextY=YM-101325EndFunction执行结果:N=2A0=353H=5E=0.1IXABC10.50820.79362788.51-52.3620.49220.90653096.52-53.67No.TK1365.5002364.2503364.8754365.1885365.0316364.953源程序中的主要

16、变量:N混合溶液的组分数目A0初值,单位KX一维数组,液相组分的摩尔分数A、B、C一维数组,安托因方程中的物质常数,单位分别为Pa、K、KT混合溶液的正常沸点,单位KX1、X2、X0、Y1、Y2、Y0、E、H等与图2-5相同。源程序前半部分采用了逐步扫描法确定根区间端点,实际上,由于理想溶液的沸点必位于两纯组分的沸点之间,故可选用两纯组分的沸点作为二分法的迭代初值。只需在该例源程序中添加上每组分的安托因常数数据(Array函数),即可用于多组分混合理想溶液的沸点计算。2.2.2迭代法2.2.2.1方法概述迭代法是一种重要的逐次逼近方法。这种方法用某个固定公式反复校正根的近似值,使之逐步精确化,

17、最后得到满足精度要求的结果。可以采用多种方法将方程式例如方程式可以通过下列不同方法变为不同形式的迭代格式:(1)等号两侧各加上x,得:(2)经移项变换后,得:(3)将方程式重新安排为这些方程式的右端都含有未知的x。如果给出根的某个近似值,并将它代入式(2-11)的右端,例如以x0为第一个近似根,则有再以x1作为第二个近似根,则依次类推,由式(2-11)可得:于是,从给定的初始近似值x0出发,按式(2-16)可以得到一个数列:如果这个数列有极限,则迭代格式(2-16)是收敛的。这时数列xk的极限就是方程式(2-11)的根。上述求非线性代数方程式数值解的方法称为迭代法。这个方法虽然简单,但根本问题

18、在于当k时,xk能否收敛于x*,也就是必须找出g(x)收敛的充分条件。这样,对于不符合收敛条件的g(x),则不必用迭代法浪费时间了。下面研究迭代格式(2-16)的收敛条件。设x*为方程x=g(x)的根,则按微分中值定理,如果f(x)的定义在闭区间a,b上而且是连续的,在开区间(a,b)内存在有限导数f(x)。那么在a与b之间至少存在一点c,满足故式中,是x*与xk之间的某一点。于是,如果存在正数q1,使得则按式(2-20)有由于x*-xk表示第k次迭代结果与根x*的误差,若令即第k+1次迭代误差Ek+1,小于或等于第k次迭代误差Ek与小于1的正数q的乘积。利用这一误差关系反复递推,最后可得考察

19、式(2-22),当k,Ek0,即g(x)收敛。结论是:如果g(x)具有连续的一阶导数,且对所有的x都存在正数q1,使得q1,则式xk+1=g(xk)对于任意初值x0均收敛。对于迭代格式xk+1=g(xk)的收敛条件,也可以用图形说明。可以把x=g(x)等号左右两侧分别表示为y=x及y=g(x)的图形(见图2-8)。由于而在y=x与y=g(x)两曲线的交点处,两曲线的函数值相等,即交点处f(x)值为零。故交点为原方程f(x)=0的解。设此交点的模坐标为x=a。当然,我们尚不知道a的准确值,但若过a的初始近似值x0对横坐标作垂线,与曲线y=g(x)相交,则其交点纵坐标A即为函数g(x0)的值。再由

20、A点对纵坐标作垂线,与线y=x相交,设其交点的横坐标为x1。因y=x是一条与横坐标成45的直线,故x1=g(x0),即x1为a的第二个近似值。按此步骤反复进行时,由于y=g(x)具有0g(x)1的特征,所以各近似值将不断逼近a,即迭代格式xk+1=g(xk)收敛。图2-9表示迭代格式具有-1g(x)0特征的迭代过程。由于满足ql,及g(x)-l时的迭代过程,这两种情况都是不能收敛的。图2-80g(x)1时的迭代情况示意图2-9-1g(x)1时的迭代情况示意 图2-11g(x)-1时的迭代情况示意由式(2-22)可知,q值越小,收敛速度越快。所以,若有多种迭代格式可供选择时,应尽可能选择在第一个

21、近似根处较小的迭代格式。例如方程式其第一个近似根x0=1.1。表2-1列出了这个方程式的5种迭代格式及其在x0=1.1处的一阶导数值。由表可知,第三种迭代格式不收敛。第五种迭代格式的收敛速度最快。表2-1不同迭代格式收敛情况序号g(x)g(x)g(1.1) 序号g(x)g(x)g(1.1)1(-3x2+4x)0.38540.46220.8855-0.18631.600下面讨论如何估计近似值xk的误差问题。按式(2-16)有式(2-24)减式(2-25)得:按微分中值定理可知:式中,为区间(xk-1,xk)中的某一点。由于迭代格式有效,故存在正数q1,且q。将此关系式代入式(2-27),得:将式

22、(2-28)代入式(2-26),得:同理如此反复递推,对任意正整数r有:利用式(2-30),可以导出计算近似值xk误差的公式。设p为任意正整数,则其中将式(2-31)中各式相加,得式(2-32)不等号右边圆括号内为一无穷递减等比级数,其各项之和可写成:故式(2-32)又可表示为令p,便有若前后两次迭代值之差等于指定的某足够小的数值:则近似值xk的误差为由于迭代格式收敛,即qEX1=(X0+1)(1/3)PrintI;:PrintFormat$(X1,#.#)OKstop=Abs(X1-X0)IfOKstopEThenX0=X1:I=I+1EndIfLoopEndSub执行结果:No.X11.3

23、572121.3308631.3258841.3249451.3247661.3247371.32472从执行结果分析,第六次迭代与第七次迭代结果之差已小于0.00001。2.2.3迭代过程的加速2.2.3.1方法概述一个迭代法的计算程序如果收敛缓慢,迭代的次数就多,计算量就大。因此,很有必要讨论加速收敛过程图2-13迭代过程的加速的方法。如图2-13所示,迭代初值为x0。第一次迭代结果为x1,x1与x*的误差用E1表示。可以设想,如果第一次迭代结果为而不是x1,便可以使迭代过程大为加速。然而,由于x*未知,故E1的数值不能确定。我们的任务是设法寻找E1的估计值,用它代替E1,以加速迭代过程。

24、假定第k+1次迭代结果为xk+1,则有按式(2-20)有若g(x)在求根范围内变化不大,近似地取为某定值q,且EX1=(X0+1)(1/3)X1=X1+(X1-X0)*Q/(1-Q)PrintTab(0);I;:PrintFormat$(X1,#.#)OKstop=Abs(X1-X0)IfOKstopEThenX0=X1:I=I+1EndIfLoopEndSub执行结果:No.X11.3224021.3247431.3247241.32472由执行结果可知,第三次迭代与第四次迭代结果之差已小于0.00001。可见达到同一准确度时,用迭代-加速格式的迭代次数大大小于简单迭代法。2.2.4牛顿法牛

25、顿法亦称牛顿-拉福森法。由于这个方法的计算结果颇佳,而计算过程亦较简单,故被普遍地使用。2.2.4.1方法概述牛顿法的核心内容是通过泰勒级数将非线性方程式转化为线性方程式,然后用迭代法求解。设方程式f(x)=0的近似根为x0,则f(x)对x=x0的泰勒级数展开式为若x=x0+h更接近于准确根,代入上式得:忽略h高于一次的幂且令f(x0+h)=0,可得:把上式中的x作为原方程的一个新的近似根,则有如下迭代格式:式(2-44)称为牛顿公式。下面用图形对牛顿公式作进一步说明。图2-15中的曲线方程为y=f(x),过曲线上k点作切线,与横坐标相交于xk+1。此切线的斜率与截距分别为故切线方程式为切线与

26、横坐标交点为方程式的根,亦即牛顿公式中新的近似根。因为牛顿公式是一种迭代格式:故有必要研究其收敛性。g(x)的一阶导数为设x*为f(x)=0的一个单根,便有符合迭代格式的收敛性条件:所以在单根x*附近,牛顿公式恒收敛,而且收敛速度很快。但是必须注意,起始值若不在根的附近,牛顿公式不一定收敛。所以,在实际使用中,牛顿法最好与逐步扫描法结合起来,先通过逐步扫描法求出根的近似值,然后用牛顿公式求其精确值,以发挥牛顿法收敛速度快的优点。此外,牛顿法的准确度与根附近的斜率有关。若根附近的斜率f(x)很小,则在计算f(x)/f(x)时,由于误差被放大而无法达到很高的准确度。计算过程要用到的一阶导数f(x)

27、,有时难于求得,这也是牛顿法的一个缺陷。使用牛顿法时,通常仍用条件来控制迭代过程结束,式中为指定的足够小的数值。2.2.4.2程序框图图2-16是先用逐步扫描法求根的近似值,再用牛顿法求精确值的通用计算程序框图。图2-15牛顿法图2-16牛顿法通用计算程序框图程序框图中的主要变量:A初值H步长E前后两次迭代结果的允许误差X0迈步起始值及迭代起始值X1迈步终值及迭代终值Y0与X0对应的函数值Y1与X1与应的函数值程序的前半部分,是用逐步扫描法求根的近似值;后半部分,是用牛顿法求根的精确值。2.2.4.3计算实例例2-7用牛顿法求解例2-5,并与迭代法求得的结果比较。解已知牛顿公式为xk+1=xk

28、-对本题由例2-5可知,原方程式在1.5附近有一个根,故以1.5为迭代初值。源程序:*Example2-7-Eg2-7.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()ClsX0=1.5:E=0.00001PrintTab(0);No.;Tab(9);X:I=1OKstop=1DoWhileOKstopEX1=X0-(X0*X0-1)*X0-1)/(3*X0*X0-1)PrintTab(0);I;:PrintFormat$(X1,#.#)OKstop=Abs(X1-X0)IfOKstopEThenX0=X1:I=I+1EndIfLoopEndSub执行结果

29、:No.X11.3478321.3252031.3247241.32472与迭代法比较,可以看出牛顿法的收敛速度甚快。例2-8CO与H2的混合气体在活性氧化锌为催化剂的情况下进行反应以生产CH3OH:CO+2H2CH3OH已知平衡常数Kp(1/Pa)Kp与温度T/K的关系式为试计算温度为365、压力为300105Pa、原料气体中CO与H2的比分别为12和32时,平衡气体中各组分的体积分数。解合成甲醇的反应是可逆、放热和摩尔数减少的反应。由于反应的可逆性,通常产品产率较低。首先要求出反应达到平衡时生成CH3OH的物质的量x。起始时物质的量0平衡时物质的量-x-2xx若令起始反应物的总物质的量为1

30、,即+=1则反应达平衡时体系的总物质的量为-x+-2x+x=1-2x故平衡时体系中各物质的摩尔分数yi为各物质的分压pi为式中,P为体系总压。将式(2-48)代入式(2-47)得式(2-49)经整理后得:令则式(2-50)写成:用D表示原料气中CO与H2的体积比,则解式(2-52)得:用牛顿法求解方程式(2-51),其迭代格式为选择初值A=0,步长H=0.03,用逐步扫描法先求出根的近似值,再令允许误差E=0.001,用牛顿法求根的精确值。源程序:*Example2-8-Eg2-8.frm*DefDblA-H,K-ZPrivateSubCommand1_Click()ClsD=InputBox

31、(D=)PrintTab(5);*DATA*A=0:H=0.03:E=0.001:T=365:P=30000000PrintT=;Format$(T,#);:Print(C)PrintP=;Format$(P/1000,#);:Print(kPa)PrintD=;Format$(D,#.#)PrintTab(4);*Result*T=T+273.2KP=3925/T-9.84*Log(T)/Log(10)-0.0034*T+9.8KP=10(KP)NC=D/(1+D):NH=1/(1+D)S=(1+KP*P*P*NH*(4*NC+NH)/(4*(1+KP*P*P)Q=KP*P*P*NC*NH*

32、NH/(4*(1+KP*P*P)X0=A:Y0=(X0-1)*X0+S)*X0-Q-OKstop=1DoWhileOKstop0X1=X0+H:Y1=(X1-1)*X1+S)*X1-QOKstop=Y0*Y1IfOKstop0ThenX0=X1:Y0=Y1EndIfLoop-OKstop=100*EDoWhileOKstopEX1=X0-(X0-1)*X0+S)*X0-Q)/(3*X0-2)*X0+S)OKstop=Abs(X1-X0)IfOKstopEThenX0=X1EndIfLoop-YC=(NC-X1)/(1-2*X1):QX=YC*100PrintY(C0)=;Format$(QX,

33、#.#);:Print(%)YH=(NH-2*X1)/(1-2*X1):QX=YH*100PrintY(H2)=;Format$(QX,#.#);:Print(%)YP=X1/(1-2*X1):QX=YP*100PrintY(CH3OH)=;Format$(QX,#.#);:Print(%)EndSub执行结果:*DATA*T=365(C)P=30000(kPa)D=0.50*Result*Y(C0)=19.62(%)Y(H2)=39.25(%)Y(CH3OH)=41.13(%)执行结果:*DATA*T=365(C)P=30000(kPa)D=1.50*Result*Y(C0)=64.11(%

34、)Y(H2)=15.35(%)Y(CH3OH)=20.54(%)源程序中的主要变量:T温度/SP压力/PaQKP平衡常数Kp/(1/Pa)YCyCONCYHNHYPD原料气中CO与H2的体积比A、H、E、X0、X1、Y0、Y1等与框图2-16相同。例2-9已知甲苯胺的饱和蒸气压计算公式为lgP=23.8296-3480.3/T-5.081lgT其中P为压力,kPa;T为温度,K。试求甲苯胺的正常沸点。解按题意,将P=101.325kPa代入上式得方程:图2-17函数f(T)的图形则由函数f(T)的几何图形(见图2-17)可见,上述方程在(0,)内有两个根,其中一个根位于300,500区间内,另

35、一个根则位于15000,20000区间内,显然后者毫无物理意义。故取初值为T0=300K,步长取H=20K,允许误差取E=0.1,按图2-16编写源程序,其中函数及其导数计算分别采用函数子程序Y和YD,以增强程序的通用性。源程序:*Example2-9-Eg2-9.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()ClsA0=300:H=20:E=0.1PrintA0=;A0:PrintH=;H:PrintE=;EX0=A0:T=X0:Y0=Y(T)-OKstop=1DoWhileOKstop0X1=X0+H:T=X1:Y1=Y(T)OKstop=Y1*Y0

36、IfOKstop0ThenX0=X1:Y0=Y1EndIfLoop-PrintNo.TKOKstop=100*E:I=0:T=X0:YD0=YD(T)DoWhileOKstopEI=I+1X1=X0-Y0/YD0T=X1:Y1=Y(T):YD1=YD(T)PrintFormat$(I,#);PrintFormat$(X1,#.#)OKstop=Abs(X1-X0)X0=X1:Y0=Y1:YD0=YD1Loop-EndSub*FunctionY(T)*Y=21.8239-3480.3/T-5.081*Log(T)/Log(10)EndFunction*FunctionYD(T)*YD=3480.

37、3/T/T-2.2067/TEndFunction执行结果:A0=300H=20E=0.1No.TK1406.0312406.1403406.140源程序中变量YD0、YD1对应与X0、X1的导数值,其他变量同框图2-16。由计算结果可见,第二次迭代结果就满足要求。此外,如取初值为T0=15000K、步长取H=50K、允许误差取E=0.1,仍可得到另一解:A0=15000H=50E=0.1No.TK118085.178218085.210显然此解毫无物理意义。可见,进行工程计算时,必须事先寻找根所在的区间,以选择合理的初值。事后还必须对所得计算结果进行分析,这样才能获得合乎实际的合理性数据。2

38、.2.5弦截法牛顿迭代法收敛速度快,但它要求计算f(x)的值。在科学与工程计算中,常会碰到f(x)不易计算或算式复杂而不便计算的情况。下面介绍一种既具有牛顿法恒收敛和收敛速度快的优点,又不用求导数f(x)的弦截法。2.2.5.1方法概述弦截法的基本思想与牛顿法相似,即将非线性函数f(x)线性化后求解。两者的差别在于弦截法实现函数线性化的手段采用的是两点间的弦线,而不是某点的切线。设方程式f(x)=0的一个近似根为x0,另一个近似根为x1。若用过点x0,f(x0)和x1,f(x1)间的弦线代替曲线y=f(x)(见图2-18),得直线方程为令式(2-54)中f(x)=0可得将上式中的x作为方程f(

39、x)=0的又一新的近似根,从而有如下迭代格式:即弦截法的迭代计算公式。显而易见,上式即为用差商代替牛顿公式(2-44)中的导数所得的结果。弦截法的几何意义就是通过点Pk-1和Pk的弦割线与x轴的交点xk+1不断逼近真实根x*的过程,如图2-18所示。与牛顿法只需给出一个初值不同,弦截法需要给出两个迭代初值x0和x1。如果与逐步扫描法结合起来,则最后搜索的区间的两个端点值常可作为初值x0和x1。弦截法虽比牛顿法收敛速度稍慢,但在每次迭代中只需计算一次函数值,又不必求函数的导数,且对初值x0和x1要求不甚苛刻,是工程计算中常用的有效算法之一。2.2.5.2程序框图图2-19是先用逐步扫描法求根的近

40、似值,再用弦截法求精确值的通用计算程序框图。图2-18弦截法的几何意义图2-19弦截法通用计算程序框图程序框图中的主要变量:A初值H步长E前后两次迭代结果的允许误差X0迈步起始值及迭代起始值X1迈步终值及迭代起始值Y0与X0对应的函数值Y1与X1与应的函数值X2迭代终值框图的前半部分,是用逐步扫描法求根的近似值;后半部分,是用弦截法求根的精确值。2.2.5.3计算实例例2-10用弦截法求解例2-9,并与牛顿法求得的结果比较。解如例2-9所述,取初值为T0=300K、步长取H=20K、允许误差取E=0.1,按图2-19编写源程序,其中函数计算采用函数子程序Y。源程序:*Example2-10-Eg2-10.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()ClsA0=300:H=20:E=0.1PrintPrintA0=;A0:PrintH=;H:PrintE=;EX0=A0:T=X0:Y0=Y(T)-OKstop=1DoWhileOKstop0X1=X0+H:T=X