非线性方程求解

非线性方程求解第1页，共32页，2023年，2月20日，星期四4.1化工实际问题的提出

求解非线性方程是化工设计及模拟计算中必须解决的一个问题。与线性方程相比，非线性方程问题无论是从理论上还是从计算公式上，都要比线性方程复杂的多。对于一般的非线性f(x)=0，计算方程的根既无一定章程可循也无直接方法可言。例如，求解高次方程组7x6-x3+x-1.5=0的根，求解含有指数和正弦函数的超越方程ex-sin(x)=0的零点。解非线性方程或非线性方程组也是计算方法中的一个主题。一般地，我们用符号f(x)来表示方程左端的函数，方程的一般形式表示为f(x)=0，方程的解称为方程的根或函数的零点。通常，非线性方程的根不止一个，而任何一种方法只能算出一个根。因此，在求解非线性方程时，要给定初始值或求解范围。而对于具体的化工问题，初值和求解范围常常可根据具体的化工知识来决定。常见的雷诺数和摩擦系数关系方程在雷诺数低于4000时有以下关系式：

(4-1)4.14.84.74.54.34.2总目录4.94.64.4第2页，共32页，2023年，2月20日，星期四4.1化工实际问题的提出这是一个典型的非线性方程。我们在管路设计中经常碰到。当我们已知雷诺数Re，如何根据公式（4-1）求出摩擦系数λ，这是我们在管路设计中必须首先解决的问题。对于方程（4-1）而言，无法用解析的方法求出摩擦系数，只能用数值求解的方法。如用在下面即将介绍的松弛迭代法，假设：

则利用松弛迭代公式可得：经11次迭代可得摩擦系数为0.07593。同样，在n个组分的等温闪蒸计算中，通过物料和相平衡计算，我们可得到如下非线性方程：

4.14.84.74.54.34.2总目录4.94.64.4第3页，共32页，2023年，2月20日，星期四4.1化工实际问题的提出在方程（4-3）中只有α是未知数，ki为相平衡常数，zi为进料组分的摩尔浓度，均为已知数。和上面的情况一样，方程（4-3）也无法直接解析求解，必须利用数值的方法，借助于计算机方可精确的计算。对于这个问题的求解，可利用我们下面介绍的牛顿迭代法进行计算，也可利用其他迭代公式进行计算，如采用牛顿迭代公式，则可以得到如下的具体迭代公式：

（4-4）

饱和蒸气压是我们经常要用到的数据，虽然我们可以通过实验测量来获取饱和蒸气压的数据，但我们通常利用前人已经测量得到的数据或回归的公式来获取，这可以减轻我们大量的基础实验工作。公式（4-5）是一种常用的饱和蒸气压计算公式：(4-5)

其中p为饱和蒸气压，单位为mmHg，T为温度，单位为K，A、B、C、D为已知系数。要想得到某一温度下的饱和蒸气压，直接利用公式（4-5）是无法得到的。因为公式（4-5）两边都有未知变量，并且无法用解析的方法求解，必须用数值计算的方法求解。通过上面的一些例子，我们可以发现，如果没有适当的手段和办法来求解非线性方程，那么化学化工中的许多研究、设计等工作将无法展开，这势必影响化学化工的发展，下面我们将介绍一些实用的非线性方程求解方法，并提供计算机程序。

4.14.84.74.54.34.2总目录4.94.64.4第4页，共32页，2023年，2月20日，星期四4.2实根的对分法

4.2.1使用对分法的条件

4.2.2对分法求根算法

4.2.3对分法VB程序清单

4.14.84.74.54.34.2总目录4.94.64.4第5页，共32页，2023年，2月20日，星期四4.2.1使用对分法的条件对分法或称二分法是求方程近似解的一种简单直观的方法。设函数f(x)在[a,b]上连续，且f(a)f(b)<0，则f(x)在[a,b]上至少有一零点，这是微积分中的介值定理，也是使用对分法的前提条件。计算中通过对分区间，逐步缩小区间范围的步骤搜索零点的位置。如果我们所要求解的方程从物理意义上来讲确实存在实根，但又不满足f(a)f(b)<0，这时，我们必须通过改变a和b的值来满足二分法的应用条件。

4.14.84.74.54.34.2总目录4.94.64.4第6页，共32页，2023年，2月20日，星期四4.2.2对分法求根算法

计算f(x)=0的一般计算步骤如下：

1、输入求根区间[a,b]和误差控制量ε，定义函数f(x)。

2、判断:如果f(a)f(b)<0则转下，否则，重新输入a和b的值。

3、计算中点x=(a+b)/2以及f(x)的值分情况处理（1）|f(x)|<ε：停止计算x*=x，转向步骤4（2）f(a)f(x)<0：修正区间[a,x]→[a,b]，重复3（3）f(x)f(b)<0：修正区间[x,b]→[a,b]，重复34、输出近似根x*。右图给出对分法的示意图。

x3=(x0+x2)/2

x2=

(x0+x1)/2

x0

x3

x1

x1

4.14.84.74.54.34.2总目录4.94.64.4第7页，共32页，2023年，2月20日，星期四4.2.3对分法VB程序清单

PrivateSubCommand1_Click()Dimx1,x2,x,y1,y2,y,eer80x1=InputBox("x1")x2=InputBox("x2")eer=inputbox(“eer”)y1=f(x1)y2=f(x2)Ify1*y2<0ThenGoTo100ElsePrint"pleaserepeatinputx1andx2"GoTo80EndIf100x=(x1+x2)/2y=f(x)IfAbs(y)<=0.001ThenPrint"thefunctionrootis";xPrint"y=";yElseIfy1*y<0Thenx2=xy2=yGoTo100Elsex1=xy1=yGoTo100EndIfEndIfEndSubPublicFunctionf(x)Dimyy=x^3+x^2-1f=yEndFunction4.14.84.74.54.34.2总目录4.94.64.4第8页，共32页，2023年，2月20日，星期四4.2.3对分法求解实例用对分法求在区间[1,2]之间的根。解：(1)f(1)=-2.8，f(2)=0.3，由介值定理可得有根区间[a,b]=[1,2]。

(2)计算x2=(1+2)/2=1.5，f(1.5)=-0.45，有根区间[a,b]=[1.5,2]。

(3)计算x3=(1.5+2)/2=1.75，f(1.75)=0.078125，有根区间[a,b]=[1.5,1.75]。一直做到|f(xn)|<ε(计算前给定的精度)或|a-b|<ε时停止。详细计算结果见表4-1。对分法的算法简单，然而，若f(x)在[a,b]上有几个零点时，如不作特殊处理只能算出其中一个零点；另一方面，即使f(x)在[a,b]上有零点，也未必有f(a)f(b)<0。这就限制了对分法的使用范围。对分法只能计算方程f(x)=0的实根。

对于多个零点的方程，我们可以通过将给定的区间[a,b]进行细分，然后在细分后的区间内用二分法分别求解，从而得到多个零点。例如求方程在0-30内的所有根。需要对二分法进行以下处理：即先给定一个a，本例中为0，然后不断增加，直到找到一个b，使f(a)f(b)<0，调用二分法，计算在[a，b]范围内的根，然后将b作为a，重复上面的工作，直到计算范围超出30为止。

4.14.84.74.54.34.2总目录4.94.64.4VB调用第9页，共32页，2023年，2月20日，星期四4.3直接迭代法

对给定的方程f(x)=0，将它转换成等价形式：。给定初值x0，由此来构造迭代序列，k=1,2,…，如果迭代收敛，即有,则就是方程f(x)=0的根。在计算中当小于给定的精度控制量时，取为方程的根。例如，代数方程x3-2x-10=0的三种等价形式及其迭代格式如下：对于方程构造的多种迭代格式，怎样判断构造的迭代格式是否收敛？收敛是否与迭代的初值有关？根据数学知识，我们可以直接利用以下收敛条件：

1、

当有

2、

在[a,b]上可导，并且存在正数L<1，使任意的，有则在[a,b]上有唯一的点满足，称为的不动点。而且迭代格式对任意初值均收敛于的不动点，并有下面误差估计式：

(4-6)

要构造满足收敛条件的等价形式一般比较困难。事实上，如果为f(x)的零点，若能构造等价形式，而，由的边疆性，一定存在的邻域，其上有，这时若初值迭代也就收敛了。由此构造收敛迭代格式有两个要素，其一，等价形式应满足；其二，初值必须取自的充分小邻域，其大小决定于函数f(x)，及做出的等价形式。4.14.84.74.54.34.2总目录4.94.64.4第10页，共32页，2023年，2月20日，星期四2.3直接迭代法例：求代数方程x3-2x-5=0，在x0=2附近的零点。

解：1）x3=2x+5构造的迭代序列收敛。取x0=2，则

准确的解是x=2.09455148150。2）将迭代格式写为

迭代格式不能保证收敛,但并不一定不收敛。VB程序界面：

4.14.84.74.54.34.2总目录4.94.64.4第11页，共32页，2023年，2月20日，星期四【例4-3】（V312）试编程计算0.01mol/L的HAc(Ka＝1.75E-5）溶液的pH值。

在新建窗体上添加一个Command1命令框（在属性窗口将Command1的属性Caption改

为计算），对命令框对象Command1的事件Click、PrivateFunctionf(X）和PrivateSubITR（XO，E，x）通用过程输入如下代码。在窗体“通用声明”段用Dim定义KA，C[能被本模块的PrivateFunctionf（x）过程存取］。

DimKA，CPrivateSubCommand1_Click）VB调用

C=0.01KA＝0.0000175XO=Sqr（K*C）

E＝1.0E-5CallITR（X0,E,x）

Print"pH=”；Int（－Log（x）／Log（10）*100＋0．5）／100EndSubPrivateSubITR（X0,E,X）

x＝XODoXO＝XX=f(XO）

LOOpWhileAbs（X—XO）＞EEndSubPrivateFunctionf(x）

f=Sqr（KA*（C－X））

EndFunction第12页，共32页，2023年，2月20日，星期四4.4松弛迭代法

有些非线性方程或方程组当用上一节中的直接迭代法求解时,迭代过程是发散的。这时可引入松弛因子,利用松弛迭代法。通过选择合适的松弛因子,就可以使迭代过程收敛。松弛法的迭代公式如下:

（4-7）由上式可知，当松弛因子ω等于1时，松弛迭代变为直接迭代。当松弛因子ω大于1时松弛法使迭代步长加大,可加速迭代，但有可能使原来收敛的迭代变成发散。当0<ω<1时,松弛法使迭代步长减小,这适合于迭代发散或振荡收敛的情况,可使振荡收敛过程加速。当ω<0时,将使迭代反方向进行,可使一些迭代发散过程收敛。松弛法是否有效的关键因子是松弛因子ω的值能否正确选定。如果ω值选用适当,能使迭代过程加速,或使原来不收敛的过程变成收敛。但如果ω值选用不合适，则效果相反,有时甚至会使原来收敛的过程变得不收敛。松弛因子的数值往往要根据经验选定,但选用较小的松弛因子,一般可以保证迭代过程的收敛。4.14.84.74.54.34.2总目录4.94.64.4第13页，共32页，2023年，2月20日，星期四4.4松弛迭代法例:用松弛迭代法求解下面非线性方程组,并分析松弛因子对迭代次数及收敛过程的影响。已知迭代初值x和y均为0,收敛精度ε=0.001

解:取以下迭代表达式:

若取松弛因子为1.1则其迭代过程如左表。若改变松弛因子,迭代过程及迭代所需的次数亦将发生变化，详见右表。由右表数据可知,当松弛因子小于1时,增大松弛因子,可加速迭代过程,减少迭代次数,但当松弛因子大于1时,迭代次数反而增加,当松弛因子达到1.56

时,迭代过程分散。

4.14.84.74.54.34.2总目录4.94.64.4第14页，共32页，2023年，2月20日，星期四此法是一种迭代加速方法。在图4-3中,从x0开始,曲线y=φ(x)和直线y=x之间的阶梯形折线为单变量的直接迭代过程。其迭代步长为每个阶梯的水平距离。若利用图4-3中曲线φ(x)上的两个点作一直线,通过求该直线和y=x的交点C的横坐标xW作为新的迭代点,这样的迭代计算过程就是韦格斯坦法。即韦格斯坦法需要在已知两点的前提下迭代计算第三点。一般第一点为人为设定,第二点利用直接迭代计算而得,则第三点就可以用韦格斯坦法迭代求解。已知A点的坐标为(x1,y1),B点的坐标为(x2,y2)。则过A、B两点的直线方程为：求该直线和y=x的交点可得:

韦格斯坦法的一般计算通式为

4.5韦格斯坦法

A

y1

y0

φ(x)Bx0x1x2x*xw图4.-3韦格斯坦法迭代示意图4.14.84.74.54.34.2总目录4.94.64.4第15页，共32页，2023年，2月20日，星期四4.5韦格斯坦法由上述公式可知,韦格斯坦法也是一种松弛法,其松弛因子为

一般情况下,当1>k>0时,迭代过程为单调收敛过程。当-1<k<0时,迭代过程为振荡收敛过程,但当k=1时,收敛将发散,故在编程计算时应注意当k=1时则取k=0进行计算。韦格斯坦法求解方程x3-2x2+5x-4=0根的QB程序见课本4.14.84.74.54.34.2总目录4.94.64.4第16页，共32页，2023年，2月20日，星期四4.6牛顿迭代法

4.6.1牛顿法的理论推导4.6.2牛顿法的几何意义4.14.84.74.54.34.2总目录4.94.64.4第17页，共32页，2023年，2月20日，星期四对方程f(x)=0可构造多种迭代格式，牛顿迭代法是借助于对函数f(x)=0的泰勒展开而得到的一种迭代格式。将f(x)=0在初始值x0做泰勒展开得：取展开式的线性部分作为的近似值，则有：设则令类似地，再将f(x)=0在x1作泰勒展开并取其线性部分得到：一直做下去得到牛顿法的迭代格式：牛顿迭代格式对应于f(x)=0的等价方程为：若b是f(x)的单根时，，则有，只要初值x0充分接近b，牛顿迭代都收敛。牛顿迭代是二阶迭代方法。可以证明，b为f(x)的a重根时，迭代也收敛，但这是一阶迭代，收敛因子为，若这时取下面迭代格式，它仍是二阶方法：

4.6.1牛顿法的理论推导

4.14.84.74.54.34.2总目录4.94.64.4第18页，共32页，2023年，2月20日，星期四以为斜率作过(x0，f(x0))点的直线，即作f(x)在x0的切线方程：令y=0，则在x1处的切线与x轴的交点x1，即：再作f(x)在x1处的切线，得交点x2，逐步逼近方程的根b。如图4-4所示。在区域[x0,x0+h]的局部“以直代曲”是处理非线性问题的常用手法。在泰勒展开中，截取函数展开的线性部分替代f(x)。

yxx2x1x04.6.2牛顿法的几何意义

图4-4牛顿切线法示意图

4.14.84.74.54.34.2总目录4.94.64.4第19页，共32页，2023年，2月20日，星期四4.6.2牛顿法的几何意义例：用牛顿迭代法求方程f(x)=x3-7.7x2+19.2x-15.3，在x0=1附近的零点。解：

计算结果列于表4-4中。

比较表4-1和表4-4的数值，可以看到牛顿迭代法的收敛速度明显快于对分法。牛顿迭代法也有局限性。在牛顿迭代法中，选取适当迭代初始值x0是求解的前提，当迭代的初始值x0在某根的附近时迭代才能收敛到这个根，有时会发生从一个根附近跳向另一个根附近的情况，尤其在导数数值很小时，如图4-5所示。x2yx1xx0x2x0x3x1如果f(x)=0没有实根，初始值x0是实数，则迭代序列不收敛。图4-6给出迭代函数f(x)=2+x2，初始值x0=2的发散的迭代序列。图4-5图4-64.14.84.74.54.34.2总目录4.94.64.4第20页，共32页，2023年，2月20日，星期四【例4-4】(V214）1mL0.2mol/LHCl溶液中含Cu2+5mg，若在室温及1.013X105Pa下通入H2S气体达饱和（0.lmol/L)则析出CuS沉淀，达平衡时，求残留在溶液中的［Cu2+］？

在新建窗体上添加一个Commandl命令框（在属性窗口将Commandl的属性Caption改为计算），对命令框对象Commandl的事件Click、PivateSubNT（X0，E，x）、PrivateFunctionf(x)和PrivateFunctiong(x)通用过程输入如下代码。在窗体“通用声明”段用Dim定义C1，C2,a,K变量为双精度型[能被本模块的PrivateFunctionf(x)和PrivateFunctiong(x)过程存取］。

VB程序第21页，共32页，2023年，2月20日，星期四

DimC1，C2，a，KAsDoublePrivateSubCommand1_Click（）

C1＝0.07874C2＝0.1a＝0.2E＝0.00001K＝11330000000000＃

CallNT（X0,E,x）

PrintxPrintFormat＄(x,”0.000e-00”）

endsubPrivateFunctionf(x)f=4*x*x-(4*(a+2*C1)+K*C2)*x+(4*C1*C1+4*a*C1+a*a)EndFunctionPrivateFunctiong(x)g=8*X-(4*(a+2*C1)+K*C2)EndFunctionPrivateSubNT（XO，E，x）

X＝XODOX0=xx=x0-f(x0)/g(x0)loopwhileAbs(x-x0)>eEndsub第22页，共32页，2023年，2月20日，星期四4.7割线法

在牛顿迭代格式中：

，用差商代导数，并给定初始值x0和x1

，那么迭代格式可写成如下形式：上式称为割线法。用割线法迭代求根，每次只需计算一次函数值，而用牛顿迭代法每次要计算一次函数值和一次导数值。但割线收敛速度稍慢于牛顿迭代法，割线法为1.618阶迭代方法。做过两点(x0,f(x0))和(x1,f(x1))的一条直线（弦），该直线与x轴交点就是生成的迭代点x2，再做过(x1,f(x1))和(x2,f(x2))的一条直线，x3是该直线与x轴的交点，继续做下去得到方程的根f(a)=0，如图4.4所示。例4.5：用割线法求方程的根，取x0=1.5，x1=4.0。解：

计算结果列于表

VB程序4.14.84.74.54.34.2总目录4.94.64.4第23页，共32页，2023年，2月20日，星期四【例4-6】(V316）计算500mL、1.0mol/LNa2S2O3溶液可以溶解多少克AgBr。在新建窗体上添加两个标签框Label1、Label2（在属性窗口将Label1的属性Caption改为割线法，属性Autosize设为True；在属性窗口将Label2的属性Caption置为空）和一个Command1命令框（在属性窗口将Command1的属性Caption改为计算），对命令框对象Command1的事件Click、PrivateSubGX（XO，H，E，XZ）和PrivateFunctionf(x)通用过程输入如下代码，在窗体“通用声明”段用Dim定义c，k变量为单精度型[能被本模块的PrivateFunctionf（x）过程存取］。VB程序第24页，共32页，2023年，2月20日，星期四Dimc,kassinglePrivatesubcommand1_click()H=0.01C=1#K1=29000000000000#Ksp=0.00000000000077K=k1*kspX0=val(inputbox(“x0=“))Callgx(x0,h,e,x)X=int(x*100+0.5)/100S=500*x/1000*188W$=chr(13)+chr(10)Label2.caption=“x=“&x&w$&”s=“&s&”g”PrintEndsubPrivatefunctionf(x)F=(4*k-1)*x*x-4*k*c*x+k*c*cEndfunction第25页，共32页，2023年，2月20日，星期四PrivateSubgx(x0,h,e,x2)bx0=f(x0)Ifbx0>0ThenX1=x0-hElseX1=x0+hEndIfDox2=X1-f(X1)/(f(X1)-f(x0))*(X1-x0)IfAbs(x2-X1)/x2<eThenExitDox0=X1X1=x2LoopEndSubnext第26页，共32页，2023年，2月20日，星期四4.8非线性方程组的牛顿方法

为了叙述的简单，我们以解二阶非线性方程组为例演示解题的方法和步骤，类似地可以得到解更高阶非线性方程组的方法和步骤。设二阶方程组

其中x,y为自变量。为了方便起见，将方程组写成向量形式：将在（x0,y0）附近进行二元泰勒展开，并取其线性部分，得到下面方程组：令则有

4.14.84.74.54.34.2总目录4.94.64.4第27页，共32页，2023年，2月20日，星期四4.8非线性方程组的牛顿方法如果再将原方程组在u1处进行二元泰勒展开，并取其线性部分，得到下面方程组：解出得出

继续做下去，每一次迭代都是一个方程组

4.14.84.74.54.34.2总目录4.94.64.4即为止。第28页，共32页，2023年，2月20日，星期四4.8非线性方程组的牛顿方法例4.7：求解下面非线性方程组取初始值解：解方