非线性方程组求解

非线性方程(组)求解2022年12月5日星期一2引言n在945.36kPa（9.33atm）、300.2K时，容器中充以2mol氮气，试求容器体积。已知此状态下氮气的P-V-T关系符合范德华方程，其范德华常数为a4.17atmL/mol2，b0.0371L/moln数学模型：范德华方程变形可得，这是关于V的非线性方程 22()()()0anf VpVnbnRTV2022年12月5日星期一3本章内容n2.1非线性方程（组）数值求解基本原理n2.2 Matlab函数n2.3 多项式求根函数rootsn2.4 非线性方程求解函数，fzero和fsolven2.5 关系和逻辑运算n2.6 Matlab程序流程控制2022年12月5日星期一4非线性方程n与线性方程相比，非线性方程问题无论从理论上还是从计算公式上都要复杂得多 n非线性方程一般形式n方程的解称为方程的根或函数的零点 n对于高次代数方程，当次数4时，则没有通解公式可用，对于超越方程既不知有几个根，也没有同样的求解方式。实际上，对于n3代数方程以及超越方程都采用数值方法求近似根 0 xf2022年12月5日星期一5数值法求根n首选要给出一个初始猜测解，然后通过各种迭代格式使其逐次逼近准确解n初值好坏对迭代收敛性有很大影响，因此初值的选取很重要n对于有专业背景的问题，初值可以按条件选择，对于没有经验的问题，可以用图解法和计算机试算搜索法初估近似解 2022年12月5日星期一62.1.1 逐步扫描法求根的初始近似逐步扫描法求根的初始近似值值n用数值法求方程的根可分为两步，首先要用数值法求方程的根可分为两步，首先要找出根的某个近似值，又称为找出根的某个近似值，又称为“初始值初始值”，然后再采用特定算法将初始值逐步接近真然后再采用特定算法将初始值逐步接近真实值，直到获得满足要求的结果实值，直到获得满足要求的结果n逐步扫描法2022年12月5日星期一7逐步扫描法0yxaby=f(x)0yxaby=f(x)a+ha+2h2022年12月5日星期一82.1.2 求方程根的精确解求方程根的精确解n非线性方程(组)的求解一般采用迭代法进行。迭代法迭代法是一种重要的逐次逼近方法。这种方法用某个固定公式反复校正根的近似值，使之逐步精确化，最后得到满足精度要求的结果n常见的迭代算法有不动点迭代、二分法、牛顿法、弦截法、威格斯坦法（Wegstein）、抛物线法等 2022年12月5日星期一9不动点迭代法不动点迭代法 0 xf xgx 我们可以通过多种方法将方程式转化为2022年12月5日星期一10迭代方程0,02ccxcxxx2xcx xcxxcxxx2122例如方程可以转化为以下不同形式(1)(2)(3)2022年12月5日星期一11方程的根n从给定的初值x0，按上式可以得到一个数列：x0,x1,x2,xk,n如果这个数列有极限，则迭代格式是收敛的。这时数列xk的极限 就是方程的根*limkkxx2022年12月5日星期一12直接迭代法n上述求非线性代数方程式数值解的方法称为直接迭代法（或称为不动点迭代法）。这个方法虽然简单，但根本问题在于当k-时，xk是否收敛于x*，也就是必须找出收敛的充分条件2022年12月5日星期一13不动点n定义定义：函数g(x)的一个不动点（fixed point）是指一个实数P，满足P=g(P)n从图形角度分析，函数y=g(x)的不动点是y=g(x)和y=x的交点 2022年12月5日星期一14不动点迭代的图形解释不动点迭代的图形解释 p0p1p2OPPxy(p0,g(p0)(p1,p1)y=xy=g(x)p0p1p2OPPxy(p0,g(p0)(p1,p1)y=xy=g(x)p0p1p2OPPxy(p0,g(p0)(p1,p1)y=xy=g(x)p0p1p2OPPxy(p0,g(p0)(p1,p1)y=xy=g(x)2022年12月5日星期一15不动点定理不动点定理 n设有(i)g,g Ca,b,(ii)K是一个正常数，(iii)p0(a,b)，(iv)对所有x a,b，有g(x)a,bn如果对于所有x a,b,有|g(x)|K1，则迭代pn=g(pn-1)将不会收敛到P点。在这种情况下，P称为排斥（repelling）不动点，而且迭代显示出局部发散性 2022年12月5日星期一162.1.2.2 二分法二分法0yxx1x2y=f(x)x*x02022年12月5日星期一172.1.2.3 牛顿法牛顿法n牛顿法也称为牛顿-拉普森法或者切线法。由于这个方法的计算结果颇佳，而计算过程也比较简单，所以被普遍采用n牛顿法的核心内容是通过泰勒级数将非线性方程式转化为线性方程式，然后用迭代法求解 2022年12月5日星期一18牛顿公式 n设方程式 的近似根为n则 对 的泰勒级数展开式为1kkkkfxxxfx 0f x 0 x f x0 xx 2000002!xxf xf xxxfxfx 0000 xfxxxfxf000f xxxfx2022年12月5日星期一19牛顿法示意图 xfy 00,xfx 11,xfx0 x1x*xYOX1kkkkfxxxfx yyx2022年12月5日星期一20牛顿法注意事项n在单根附近，牛顿公式恒收敛，而且收敛速度很快。但是需要注意如果初始值不在根的附近，牛顿公式不一定收敛n在实际使用中，牛顿法最好与逐步扫描法结合起来，先通过逐步扫描法求出根的近似值，然后用牛顿公式求其精确值，以发挥牛顿法收敛速度快的优点 2022年12月5日星期一212.1.2.4 弦截法弦截法n牛顿迭代法收敛速度快，但它要求计算函数导数的值。在科学与工程计算中，常会碰到函数导数不易计算或者算式复杂而不便计算的情况 n弦截法的基本思想与牛顿法相似，即将非线性函数线性化后求解。两者的差别在于弦截法实现函数线性化的手段采用的是两点间的弦线，而不是某点的切线 111kkkkkkkfxxxxxfxfx2022年12月5日星期一22弦截法弦截法示意图 11kkkkkxxxfxfxf 00,xfx 11,xfx0 x1x*xYOX2x 22,xfx111kkkkkkkfxxxxxfxfx kkkkxfxfxx12022年12月5日星期一23弦截法注意事项n与牛顿法只需给出一个初值不同，弦截法需要给出两个迭代初值。如果与逐步扫描法结合起来，则最后搜索的区间的两个端点值常可作为初值n弦截法虽比牛顿法收敛速度稍慢，但在每次迭代中只需计算一次函数值，又不必求函数的导数，且对初值要求不甚苛刻，是工程计算中常用的有效计算方法之一 2022年12月5日星期一242.1.2.5 松弛迭代法松弛迭代法n有些非线性方程用前面的不动点迭代法求解时，迭代过程是发散的。这时可以引入松弛因子，利用松弛迭代法。通过选择合适的松弛因子，就可以使迭代过程收敛 1nnnnxxxx 2022年12月5日星期一25松弛法注意事项n由上式可知，当松弛因子=1时，松弛迭代法变为不动点迭代法；当松弛因子1时，松弛法使迭代步长加大，可加速迭代，但有可能使原理收敛的迭代变为发散；当01时，松弛法使迭代步长减小，这适合于迭代发散或振荡收敛的情况，可使振荡收敛过程加速；当0时，将使迭代反方向进行，可使一些迭代发散过程收敛 2022年12月5日星期一26松弛法注意事项n松弛迭代法是否有效的关键因素是松弛因子的值能否正确选定。如果值选用适当，能使迭代过程加速，或者使原来不收敛的过程变成收敛；但如果值选用不合适，则效果相反，有时甚至会使原来收敛的过程变得不收敛。松弛因子的数值往往要根据经验选定，但选用较小的松弛因子，一般可以保证迭代过程的收敛 2022年12月5日星期一272.1.2.6 威格斯坦法威格斯坦法n威格斯坦法在化工流程模拟中得到了广泛威格斯坦法在化工流程模拟中得到了广泛应用应用 n威格斯坦法是一种迭代加速方法 2022年12月5日星期一28Wegstein法的示意图和迭代公式)()()()()(11111nnnnnnnnnxxxxxxxxx2022年12月5日星期一29Wegstein法注意事项n应注意，如果x1和x2两点选择不当，则连线的斜率等于1，与直线y=x无交点，从而迭代无法进行，这就是Wegstein法应当避免的陷井。引入一个量CSC1111)()(nnnnxxxxS)()1(1nnnxCxCx2022年12月5日星期一30Wegstein法注意事项n令q1-Cn当q0时，Wegstein法退化为简单的不动点迭代n当0q0时，迭代能稳定收敛，但收敛较慢n当q0可以加速收敛，但易导致不稳定n为了加速收敛又避免不稳定，常取-5qc=1-1 0-1;r=roots(c)r=1.4656 -0.2328+0.7926i -0.2328-0.7926ipolyval(c,r(1)ans=-2.5535e-015 2022年12月5日星期一642.4 非线性方程求解函数nfzeroq一般的单个超越方程，可以采用fzero函数求解。它结合使用二分法、割线法和可逆二次内插法 nfsolveq与fzero函数只能求解单个方程的根不同，fsolve函数可求解非线性方程组的解。其算法采用的是最小二乘法 2022年12月5日星期一652.4.1 函数函数fzeronx,fval,exitflag,output=fzero(fun,x0,options,p1,p2,.)此函数的作用求函数fun在x0附件的零值点，x0是标量qfval函数在解x处的值qexitflag 程序结束情况:0，程序收敛于解；0，程序没有收敛；0，计算达到了最大次数qoutput 是一个结构体，提供程序运行的信息；output.iterations，迭代次数；output.functions，函数fun的计算次数；output.algorithm，使用的算法qoptions 选项，可用optimset函数设定选项的新值qfun可以是函数句柄或匿名函数。2022年12月5日星期一66例题例题7：n计算以下方程的根q1)求sinx在3附近的零点;q2)求cosx在1,2范围内的零点；q3)q4)n本例较简单，可直接在命令窗口输入命令求解：q1)fzero(sin,3)q2)fzero(cos,1,2)q3)fzero(x)x3-2*sin(x),1)q4)fzero(x)x3-2*x-5,1);roots(1 0-2-5)0sin23xx0523 xx2022年12月5日星期一67说明说明1)第2小题中，如果所给区间两端方程不异号，则程序出错2)除了采用匿名函数外，当然可以采用句柄函数定义函数，例如第三小题可以采用如下程序：function Cha2demo1 x=fzero(fun,1)function y=fun(x)y=x3-2*sin(x);3)初值的选择对于解有影响，不同的初值可能获得不同的解q可以根据感兴趣的解的区间确定初值范围q可以作出函数在一定范围内的曲线，直观的确定解的大致范围4)fzero不能获得多项式的多重根，尤其是复数根。而roots函数求解，则可获得所有根2022年12月5日星期一68例题例题8n求 的零点，以t为自变量，取值范围为-10t(大于)，=(大于等于)，=a2的运算结果 a1=1 2 3;4 5 6;7 8 9;a2=5;a3=a1=a2a3=0 0 0 0 1 1 1 1 1 1234567892022年12月5日星期一762.5.2 逻辑运算逻辑运算2.5.2 逻辑运算逻辑运算n两个逻辑量之间可以进行“与”、“或”和“非”三种基本逻辑运算及由它们组合而成的其它运算nMatlab中的逻辑操作符包括：&(逻辑与)，|(逻辑或)，(逻辑非)，xor（异或），&（先决与）和|（先决否）n逻辑运算的结果仍然是逻辑量0（假）或1（真）n维数相同的矩阵进行逻辑运算时，定义为它们对应运算逻辑量间的逻辑运算结果nMatlab允许数与矩阵间进行逻辑运算，规则与关系运算相同，是数与矩阵各个运算间的逻辑运算2022年12月5日星期一77逻辑运算真值表 逻辑量及其运算真值（逻辑量）A1100B1010A0011B0101A&B1000A|B1110 xor(A,B)0110A&B1000A|B11102022年12月5日星期一78例题例题13：求数值矩阵的逻辑求数值矩阵的逻辑“非非”，a1和和0的的“异或异或”解：解：键入a1=1 0-5 0;3-2 0 6;0 0 5 7;9 2 1 9;a2=a1,a3=xor(a1,0)回车得到结果a2=0 1 0 1 0 0 1 0 1 1 0 0 0 0 0 0a3=1 0 1 0 1 1 0 1 0 0 1 1 1 1 1 12022年12月5日星期一792.5.3 运算符的优先级运算符的优先级nMatlab表达式中可能包含多种运算符：数学运算符、关系运算符和逻辑运算符n各运算符执行的先后是根据优先级别执行的n不管运算符的位置如何，具有高优先级的运算符先执行n具有相同优先级别的运算符则按先左后右的次序执行2022年12月5日星期一80运算符的优先级优先级运算符1（）2.3+(uplus)-(uminus)4.*./*/5+-6:7=0y=sqrt(x)elseerror(Negative input not allowed)end2022年12月5日星期一862.6.2 switch多重分支结构多重分支结构 switchcaseotherwise语句的能力与ifelseend语句类似，但对多重选择的情况switch语句使代码更加易读。switch的一般格式：的一般格式：switch test_expr%测试表达式test_expr可以是标量或字符串case value statements 当test_expr值是value时，执行该语句组case value1,value2,statements 当test_expr值是value1或value2或时，执行该语句组otherwise statementsend2022年12月5日星期一87例题例题15：编程判断并显示输入自变量x的具体数值时，函数是奇数、偶数（0视为偶数，不分正负），还是其它解：解：程序如下%switchclear,x=input(x=)y=x4-19*x3+x+32;switch mod(y,2)case 1 disp(It is odd)case 0 disp(It is even)otherwise disp(It is others)end2022年12月5日星期一882.6.3 for循环结构循环结构nfor循环结构用于循环次数已知的情况，其一般格式如下：for 循环变量=表达式1（初值）：表达式2（步长）：表达式3（终值）statements（语句组）endnforend循环的执行过程是：先计算初值和终值，并把初值赋给循环变量；再判断循环变量的值是否超过了终值。若超过，则退出循环，执行end后面的语句，否则执行循环体的语句组，之后将循环变量加上一个步长，然后重复执行循环体内容，直至循环变量超过终值而退出循环为止2022年12月5日星期一89例题例题16：n编程计算1至1000的加和值解：解：编程如下clearj=0;for i=1:1000 j=j+i;endji运行后结果为j=500500i=1000n可见循环变量的最终值为表达式的终值2022年12月5日星期一902.6.4 while循环结构循环结构while cndition（表达式）statements（执行语句组）end当MATLAB执行这个whileend循环时会首先测试condition（条件表达式）。如果condition为False(零)，则直接跳出循环，执行end后面的语句。如果condition为True(非零)，则执行语句组statements，然后退回到while语句再测试条件例题例题17：求满足的最小整数m解：解：由于终止条件未知，采用while-end循环。cleara=0;n=1;while a=1e4 a=a+n2;n=n+1;endm=n-1,a2022年12月5日星期一912.6.5 trycatch语句语句ntrycatch是用于对异常进行处理的语句。把有可能引起异常的语句放在try控制块中，这样当try控制块中的statement语句引起异常时，catch控制块就可以捕获它，并针对不同的错误类型，进行不同的处理。它与C+中的trycatch语句作用一样2022年12月5日星期一922.6.6 编程中的控制语句编程中的控制语句n2.6.6.1 continue语句语句 continue语句通常置于for循环或while循环内，根据条件执行continue语句。当执行continue时，即跳出循环体中尚未执行的语句，接着进行下一次是否进行循环的判断。n2.6.6.2 break和和pause语句语句 break语句通常也置于for循环或while循环内，根据条件执行break语句，以直接退出最内层的for循环或while循环。而用pause语句将使程序暂停以等待用户按键，当用户按下任意键时程序才继续往下执行2022年12月5日星期一932.6.6.3 return语句语句n在某个函数内部执行return语句时，可立即退出该函数，并返回到调用它的函数，继续运行。当函数过程已完成每个任务并可直接返回时，return是非常有用的nreturn通常放置于函数内的一个控制结构（如if语句）内。执行该函数时，如果符合控制结构的某个条件，则调用return语句终止当前运行，并返回到调用它的函数或环境n该语句的作用与其他高级语言（如FORTRAN）的return语句的作用相同