Ch4非线性方程组及其应用

化工应用数学AppliedMathematicsforChemicalEngineering第4章非线性方程（组）及其应用化学化工学院汤吉海2概述求解方程f(x)=0的数值法的基本思想：从某个初始近似值x0出发，按照某种数值过程模式进行重复，从而逐次改进前次的结果，直到达到规定的精度，便停止这种重复，并将最后所得结果xn作为方程f(x)=0的数值解。数值法解方程必须做三件事：首先:在函数f(x)的定义域内寻找初值x0；其次:建立逐次逼近的数值过程模式；第三:规定最终近似解的精度。解的精度取决于工程计算的实际要求而与算法无关。各种算法的实质是由数值过程的模式确定的，不同的模式具有不同的计算效果和速度。初值x0则影响计算过程的速度甚至成败。3主要内容4.1非线性方程的求解二分(Bisection)法简单迭代(Iteration)法牛顿(Newton)法弦截(Secant)法抛物线法（Muller法）4.2非线性方程组的求解牛顿-拉夫森(Newton-Raphson)法4.3非线性方程（组）的MATLAB求解图解法数值求解(fzero，fsolve)4.4化工应用应用1：根据PVT关系计算摩尔体积应用2：多组分溶液体系的泡点应用3：绝热连续搅拌釜式反应器的转化率(自学)44.1.1二分(Bisection)法原理：设f(x)在区间[a,b]上连续，且有f(a)f(b)0，则由连续函数介值定理，f(x)在[a,b]内必有零点，称[a,b]为方程f(x)＝0的有根区间。二分法：取[a,b]区间的中点x0=(a+b)/2。将区间[a,b]分成两半，然后逐步搜索求解，若f(x0)与f(a)同号，则令a1=x0，b1=b;反之，若f(x0)与f(a)异号，则令a1=a，b1=x0。这样得到新的有根区间[a1,b1]，其长度[a,b]为之一半。特点：二分法计算简单，但是收敛比较缓慢；当存在多个实根时，则搜索策略要复杂很多。54.1.2简单迭代(Iteration)法迭代法是一种重要的逐次逼近方法，这种方法采用某个固定公式反复校正根的近似值，使其逐步精确化。求解方程f(x)=0的实根，将方程改写成等价形式：x=φ(x)，φ(x)是x的新函数，xk+1=φ(xk)称作迭代公式。假设x0是方程根的初始近似值，则代入迭代公式右端，求得x1，再将x1代入右端，求得x2,…，如此重复循环，即：x1=φ(x0)，x2=φ(x1)，…，xn=φ(xn-1)采用不同的方法改写求解方程则得到不同的迭代公式，求解的收敛性也不同。6简单迭代法的收敛性发散收敛74.1.3牛顿(Newton)法基本思想：将非线性方程线性化，以线性方程的解逐步逼近非线性方程的解。),2,1,0()()('1kxfxfxxkkkkNewton法是以曲线的切线与x轴的交点作为曲线与x轴的交点的近似。故Newton法又称为切线法。8从几何上看，弦截法是以曲线上两点的割线与x轴的交点作为曲线与x轴的交点的近似，故弦截法又称为割线法。4.1.4弦截(Secant)法)()()()(111kkkkkkkxxxfxfxfxx此外，割线的函数表达式恰为函数f(x)的以xk，xk-1，为节点的线性插值多项式。用弦截法需要给出两个初始值x0，x1，通常取根所在区间的端点即可。94.1.5抛物线（Muller）法如果考虑用f(x)的二次插值多项式的零点来近似f(x)的零点，就导出抛物线法。迭代公式见教材P70。此法必须已知方程f(x)=0的根的三个近似值，以这三点为节点的f(x)的二次插值多项式如图所示，抛物线法的几何意义是以过曲线上三点的抛物线与x轴的交点作为曲线与x轴交点的近似。实际计算表明，抛物线法对初值要求并不苛刻。在初值不太好的情形下常常也能收敛。它的缺点是程序较复杂。104.2非线性方程组的解法解非线性方程组的方法主要有两类：其中一类是线性化的方法，它是将方程组中的每个方程线性化得到一个线性方程组，由此构造迭代格式，求得非线性方程组的近似解。其代表是Newton法；另一类方法是将解非线性方程组问题化成优化问题，然后以最优化方法求解。最速下降法是其中最基本的方法。0),,,(0),,,(0),,,(21212211nnnnxxxfxxxfxxxf114.2牛顿-拉夫森(Newton-Raphson)法考虑一个由两个方程构成的非线性方程组：考虑方程组在估计值（x1k，x2k）处的一阶Taylor展开式：将此方程组写成如下的矩阵-向量形式：由此可以得到新的估计值：若令：则方程求解过程可以写成：牛顿－拉夫森法（Newton-Raphson）12Newton-Raphson法非线性方程组的Jacobian矩阵knnnnnkxxxfxfxfxfxfxfxF2112111')(0),,,(0),,,(0),,,(21212211nnnnxxxfxxxfxxxf计算过程的停止准则可以采用以下两种：),,,(max)(1121111knkkinikxxxfxFkikinikkxxxx111max当较大时，也可以用两次迭代的相对误差充分小作为停止准则。kx13Newton-Raphson法举例14Newton-Raphson法源程序function[x,fx]=newton(f,df,x0,TolX,MaxIter)TolFun=eps;xx1=x0;fx=feval(f,x0);fork=1:MaxIterdfdx=feval(df,xx1);%derivativefunctiondx=-dfdx\fx;xx2=xx1+dx;%IterationEqfx=feval(f,xx2);ifabs(fx)TolFun|abs(dx)TolXbreak;endxx1=xx2;endx=xx2;ifk==MaxIterfprintf('Thebestin%diterations\n',MaxIter)elsefprintf('Themaximumiterationsis%d\n',k-1)end15Newton-Raphson法举例函数定义：Jacobian矩阵定义：functiony=f42(x)y(1)=x(1)^2-2*x(1)-x(2)+0.5;y(2)=x(1)^2+4*x(2)^2-4;y=[y(1);y(2)];endfunctiondy=df42(x)dy(1,1)=2*x(1)-2;dy(1,2)=-1;dy(2,1)=2*x(1);dy(2,2)=8*x(2);endclear;x0=[2.0;0.25];TolX=1e-6;MaxIter=50;[x,err]=newton('f42','df42',x0,TolX,MaxIter)Themaximumiterationsis3x=1.90070.3112err=00164.3非线性方程（组）的MATLAB求解4.3.1代数方程的图解法采用函数ezplot(f)画出表达式f=f(x)的图，默认的范围是：-2πx2π。也可以指定范围，如ezplot(f,[min,max])画出f=f(x)的图，区域范围是：minxmax。holdoffezplot('exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5',[05])holdonline([0,5],[0,0])grid5.0)2cos(4)24sin(5.03tetettt=3.52;x=exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5x=3.1545e-00417图解法图解法由于受到空间维数的限制，对超过2维的非线性方程组的求解不太合适，需要寻求其他的求解方法。holdoffezplot('x^2+y^2-4')%画第一条曲线holdonezplot('0.75*x^3-y+0.9')%画第二条曲线grid09.075.004322yxyx184.3.2多项式型方程的求根对于多项式：可以使用roots函数求根，使用方法为：roots([anan-1an-2…a1a0])例：001111axaxaxannnn0132245xxxroots([230101])ans=-1.7246-0.3632+0.6597i-0.3632-0.6597i0.4755+0.5339i0.4755-0.5339i194.3.3非线性方程组的准解析解法例：解方程：检验结果：g=solve(eq)%最简单的调用g=solve(eq,var)%简单调用，并指定变量g=solve(eq1,eq2,...,eqn)%直接得出根g=solve(eq1,eq2,...,eqn,var1,var2,...,varn)%直接得出根，并指定变量09.075.001322yxyx[x,y]=solve('x^2+y^2-1=0','0.75*x^3-y+0.9=0')[eval('x.^2+y.^2-1'),eval('0.75*x.^3-y+0.9')]ans=[0.,-.1e-31][0.+0.*i,0.+0.*i][0.+0.*i,0.+0.*i][-.1e-31,0.][-.6e-31-.1e-30*i,-.20e-30+.3e-31*i][-.6e-31+.1e-30*i,-.20e-30-.3e-31*i]20上例结果x=.35696997189122287798839037801365.86631809883611811016789809418650+1.2153712664671427801318378544391*i-.55395176056834560077984413882735+.35471976465080793456863789934944*i-.98170264842676789676449828873194-.55395176056834560077984413882735-.35471976465080793456863789934944*i.86631809883611811016789809418650-1.2153712664671427801318378544391*iy=.93411585960628007548796029415446-1.4916064075658223174787216959259+.70588200721402267753918827138837*i.92933830226674362852985276677202+.21143822185895923615623381762210*i.19042035099187730240977756415289.92933830226674362852985276677202-.21143822185895923615623381762210*i-1.4916064075658223174787216959259-.70588200721402267753918827138837*i21符号计算示例例：解方程：3036222xayybaxx[x,y]=solve('x^2+a*x^2+6*b+3*y^2=0','y=a+x+3','x,y')x=1/2/(4+a)*(-6*a-18+2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b-3*