![薛定宇教授大讲堂(卷Ⅳ):MATLAB最优化计算](https://wfqqreader-1252317822.image.myqcloud.com/cover/152/29977152/b_29977152.jpg)
2.3 代数方程的数值求解
前面介绍的图解方法只是非线性方程组众多求解方法的一种,图解法有其明显的优势,但也有劣势。图解法只能用于求解一元或二元方程的实数根,对多元方程是不能采用图解法求解的。本节将探讨一般方程的求解思路与实用求解函数。
2.3.1 Newton–Raphson迭代方法
为简单起见,可以先探讨一元方程的求解方法。Newton–Raphson迭代方法是以英国科学家Isaac Newton与英国数学家Joseph Raphson(约1648−约1715)命名的一般方程的迭代方法。
假设一元方程是由f(x)=0描述的,且在x=x0点处函数的值f(x0)为已知的。这样,可以在(x0,f(x0))点作函数曲线的一条切线,如图2-8所示,则该切线与横轴的交点x1可以认为是找到的方程的第一个近似的根。由图2-8给出的斜率为f′(x0)的切线方程,可以得出x1的位置为
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P28_29617.jpg?sign=1738892835-yKNsVOp78qjzUdyWNop7D2YGiI9aFqJ9-0-537592d4515a6d101879031c0abba02a)
式中,f′(x0)为f(x)关于x的导函数在x0点的值。再由x1出发作切线,则可以得出x2,以x2出发找到x3,…。若已经找到了xk,则从该点出发可以搜索到下一个点。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P28_29619.jpg?sign=1738892835-fWLlx7BQnoxm9fPQhWmMKHU4R1bZl68X-0-3edba30df6bd488833cb154196af61f1)
如果|xk+1−xk|≤ε1或|f(xk+1)|≤ε2,其中,ε1与ε2为预先选定的误差容限,则可以认为xk为原方程的一个解。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P28_29620.jpg?sign=1738892835-G0dS8NptgI7lqgvHyGDhxSrs3IK07frY-0-31ebc7f6a6c3af94077bba86d1819f2e)
图2-8 Newton–Raphson迭代法示意图
定义2-7 对多元向量函数f(x),其Jacobi矩阵的定义为
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P28_29621.jpg?sign=1738892835-AT0vw6QuIiGG1M8B3C6lLT9LxgHzyGFG-0-304a5e4726fda2250faddc071a4a8af8)
定理2-3 更一般地,对多元方程f(x)=0,其中x为向量或矩阵,而非线性函数f(x)也是同维数的函数,则可以由下式迭代求出:
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P28_29624.jpg?sign=1738892835-RxT0KIUHf2FnlgkTBSgnPzgFMNIlxjo6-0-5f1cdb03e4f84148d727b8334d2b44bf)
如果||xk+1−xk||≤ε1或||f(xk+1)||<ε2,则xk为方程的根。在定义中使用了Jacobi矩阵。
根据给出的算法,编写出MATLAB的通用求解函数:
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P28_29625.jpg?sign=1738892835-IuqRO2fkZwMgn2jM67lRo8QEsiX3fq3F-0-1dda7942d874707f00ef8608d6829fd4)
该函数需要用户提供方程函数句柄f及其Jacobi矩阵的句柄d,还需要给出初值的列向量x0,并给出误差限ε。如果给出key,并令其值为1,则可以将中间搜索结果由x矩阵返回。
对单变量函数而言,如果需要提供给定函数的导函数本身就是个比较麻烦的事,使用可以考虑采用正割方法来近似函数的导数,搜索方程的根。这时,求解式(2-3-2)可以改写成
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P29_29629.jpg?sign=1738892835-Eek3EaAYqnaOPb3ugCYhzx7hyyabXn0o-0-1acaeeda3cfbc060db4a516d9726878e)
如果采用点运算,则这里的正割函数同样适用于多元方程的求解。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P29_29631.jpg?sign=1738892835-D4uAUTstIryCxbyPGLGvuqZIvdMrwlMr-0-eac5dc9f28fba7e085a52f0ca75d0aba)
例2-13 选择初值x0=10,并求出例2-10中一元超越方程的一个根。
解 先由符号运算的方式推导出给定函数的一阶导数。
>> syms x; f=exp(-0.2*x)*sin(3*x+2)+cos(x), diff(f)
可以得出函数的导数为
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P29_29633.jpg?sign=1738892835-zKdHtE7yfcbV4An4nnHPkiOdA7LDqKD6-0-686f964872d1e26d10ac896b8d844598)
有了函数及其导数,则可以用匿名函数表示它们,再设定搜索的初值为x0=10,这时调用Newton–Raphson求解算法,则可以得出方程的解搜索的中间点为[10,10.8809,11.0700,11.0593,11.0593]T,方程的解为x=11.0593,将其代入原方程,则可以得出误差为−5.1348×10−16。可以看出,利用这样的方法求解精度还是比较高的。整个求解过程的示意图如图2-9所示,可以看出经过几步迭代就可以得出方程的解。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P29_29637.jpg?sign=1738892835-J7G8hMvYrgcxHWRsei2H1wkI8VZDVGX9-0-0d143e867c191f3b7da4b777a7bac67b)
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P30_29638.jpg?sign=1738892835-NhU5mCcVzgCkTXc6Wylj3xY4D8jtdTyZ-0-e243be79aba014b00ba752c293f59ded)
图2-9 一元方程求解的中间过程
例2-14 试用正割法重新求解例2-13中方程的一个根。
解 选择两个初始点x0=9,x1=10,调用求解函数则可以得出求解过程的中间结果为x=[9,10,12.9816,11.3635,10.7250,11.0222,11.0633,11.0592,11.0593,11.0593],其中搜索过程如图2-10所示。可以看出,虽然这里给出的求解函数效率不如Newton–Raphson算法,但其优势是无须用户提供函数的导函数,所以该算法是有意义的。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P30_29641.jpg?sign=1738892835-BfgSzd8dNa3zrrSIhsyiZ4tNFVUbDm3Y-0-f5150b850bf04e32c951163e2c077bc4)
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P30_29642.jpg?sign=1738892835-Vo3sdsWxktkhOpYxHox52K5SXbzKW4HR-0-2d134aaf1a1484e78d3acc7411a6fb3a)
图2-10 一元方程求解的中间过程
例2-15 试求解例2-11中的二元方程,以(1,1)点为初值搜索到一个解。
解 由于同时含有自变量x和y,这样的方程是不能直接求解的,需要将其改写成x的方程,最简单的方法是令x1=x,x2=y,这样,原始的方程可以改写成
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P30_29646.jpg?sign=1738892835-nmuSrut7Sf422mu0sZeJbxYHodpwF5Mf-0-67be261874d8319a76f45f974e75a2b1)
Jacobi矩阵不易用手工的方法推导出来,是需要通过解析推导求解的,所以应该在符号运算框架下输入原始函数,并通过jacobian()函数计算出该矩阵。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P31_29649.jpg?sign=1738892835-TtFndbllDxjpd2xdbmO2QAAlFrKWgdrn-0-6915190fe55d02337e1507d3dbd8fe53)
可以推导出函数的Jacobi矩阵为
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P31_29651.jpg?sign=1738892835-P4SJL0U3v32hZjhYaQ4sL2HfG8sweXyH-0-47b86d927ebdced0cb1820d26d84ad99)
有了原函数与Jacobi矩阵,则可以手工写出这两个函数的匿名函数,然后调用基于Newton–Raphson算法的求解函数。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P31_29653.jpg?sign=1738892835-OI3WKGOjrHT6bioPyLdpWwKJh9mU42PC-0-c0313a5a1cef5a5ead859be1ee8b0946)
由该初值点出发得出的方程的解为x=[5.1236,−12.2632]T,中间点的个数为18。代入原方程后得出的误差矩阵范数为3.9323×10−12。从求解的结果看,确实通过该函数可以求出方程的一个根,不过从实际操作看,这样的方法未免过于复杂,由于需要推导Jacobi矩阵,并将其矩阵用匿名函数的形式手工表示出来,该过程比较容易出错。
若采用正割方法,则可以给出下面的语句,得出方程的根为x=[1.0825,−1.1737]T,代入方程得出误差为2.2841×10−14。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P31_29655.jpg?sign=1738892835-ZBNXfbvzx6J90kos2DYQIuUIHAdG1bB1-0-28ae8e5358caac603247ee5641e817b5)
虽然这时函数调用无须用户提供Jacobi矩阵,但所需的迭代步数为265,求解效率较低,所以对代数方程求解问题而言,需要更好的方法。
2.3.2 MATLAB的直接求解函数
由于上面给出的求解算法比较麻烦,需要提供的信息又比较难获取,所以在实际方程求解中应该考虑采用更好的方法。MATLAB提供了更实用的求解函数fsolve(),无须提供Jacobi矩阵的句柄,只需给出方程函数的句柄和初值,则可以直接求解任意复杂的非线性方程组,由给出的初值搜索出方程的一个根。该函数的调用格式为
x=fsolve(f,x0)
式中,f为方程函数的句柄;x0为初始向量或矩阵;f函数的维数与x0完全一致,正常情况下得出的x为方程的数值解。
该函数可以至多返回四个变量,这时完整的调用格式为
[x,F,flag,out]=fsolve(f,x0,opts)
式中,x为方程的解;F为x处方程函数的值矩阵;flag如果为正说明求解成功,out变量还将返回一些中间信息。用户还可以增加输入选项opts控制求解的算法与精度,后面将通过例子演示。
例2-16 试利用这里介绍的求解函数直接求解例2-15中的方程。
解 仍然可以使用匿名函数描述原始的方程,且无须提供Jacobi矩阵的函数句柄,直接调用求解函数即可以得出方程的解为x1=[0,2.1708]T,将解代入方程则得出误差向量的范数为3.2618×10−14。虽然这样得出的解不是例2-15中得出的解,但也是方程的一个解。此外,还可以看出,迭代步数为14次,f函数的调用次数为38,与例2-15中的结果相仿。不过该函数的优势是无须用户提供方程的导函数,使用该函数更适合于实际应用,建议采用该方法直接求解方程。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P32_29664.jpg?sign=1738892835-tmLkzs9TnDbmNReVDMkRPaq0uKx44qQh-0-50d4a9c2570648335ebe6a446bd6e7fe)
如果将初值修改为x0=[2,1]T,则搜索到的解为x1=[−0.7038,1.6617]T,该解对应的误差为f1=2.0242×10−12。
>> x0=[2; 1]; x1=fsolve(f,x0), f1=norm(f(x1))
与前面介绍的Newton–Raphson迭代法相比,求解方程的过程已经极大地简化了,由于只需描述方程本身,无须描述更复杂的Jacobi矩阵,使得方程的求解变得轻而易举,所以在实际方程求解问题中可以放心使用这样的方法。
在前面的例子中,未知自变量x与方程函数f(x)都是同维数的向量,编写出匿名函数就可以描述方程函数,就可以求解方程从而得出方程的解。如果方程f(x)=0中,x与f为同维数的矩阵,也可以直接利用fsolve()函数搜索出方程的数值解。下面以Riccati方程为例介绍矩阵型方程的求解方法。
定义2-8 Riccati代数方程的数学形式为
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P33_29670.jpg?sign=1738892835-WY8cZNroNYLr3dw4rRCUfNewGxl5Ebsz-0-fe3bb25cef3016538fae05cf61513098)
式中,每个矩阵都是n×n矩阵。
Riccati方程是以意大利数学家Jacopo Francesco Riccati(1676−1754)命名的,本意是对应于一类一阶微分方程,其原型要求B为正定矩阵,C为对称矩阵。后来因为微分方程难以求解,所以将其简化成上面的Riccati代数方程。
从数学角度看,这几个矩阵可以为任意矩阵。在控制科学领域,可以考虑采用控制系统工具箱提供的are()函数直接求取方程的数值解,不过该函数只能求出Riccati方程的一个根。如果想获得方程全部的根,则可以考虑调用vpasolve()函数直接求解。下面将给出例子演示Riccati方程与多项式矩阵方程的方法。
例2-17 试求解下面的Riccati代数方程。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P33_29675.jpg?sign=1738892835-Uk2x1hTudkMqUIkcZCdU0YryjfgAtQsS-0-d198a9e7abd9f7d28f30f73b91a93675)
解 可以先输入这几个矩阵,然后调用控制系统工具箱中的are()函数求解Riccati方程,并计算得出的误差矩阵范数。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P33_29677.jpg?sign=1738892835-oxrSsJsHIeHJipYF7UY9DTNB2yvU2HSk-0-77ffe51b8338db09a0d862418c224539)
由控制系统工具箱的are()函数可以直接得出方程的一个根如下,代入原方程可见,该解导致的误差为1.3980×10−14,说明该解的精度是比较高的。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P33_29679.jpg?sign=1738892835-CPeJEhmDs80P4Gp6GnCrYGP0SXtrn0vd-0-5a9d55982f321b856c02569a7432771c)
由于该方程是二次型方程,人们很自然就可以想到,该方程可能有多个根,如何求出其他的根呢?显然可以尝试选择一个不同的初始矩阵,例如幺矩阵,从该矩阵求解方程的根。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P33_29683.jpg?sign=1738892835-KBToZklPulgCPtpOCFZxLvY96PxFpeL9-0-d62dae399c04bd10d176f5e9a2c1441c)
得出方程的解如下,对应的误差为6.3513×10−9。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P33_29681.jpg?sign=1738892835-ulp6yBQY4kpAMTtdDjLRZKMzs1qCUm8Q-0-957beecc5d026346d56c6735d87c1ea6)
显然,这时得出的解与are()函数得出的解是不一致的。该函数返回的f1矩阵就是方程解处的函数值,即f(X1)。另外,在这个例子中返回的flag值为1,因为它为正数,所以表示求解成功。另外返回的cc信息包括如下内容,表明迭代步数为11,函数调用次数为102,说明求解效率还是很高的。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P34_29686.jpg?sign=1738892835-IfmKeVJh47oHxO2Qt4FRAXOfzsMO3S8B-0-c735f460f8dd4d2911a35152a01dd410)
2.3.3 求解精度的设置
从上面例子得出的解看,误差偏大。有没有什么办法控制误差的大小呢?
前面在定理2-3中描述了一般迭代方法收敛的条件,给出了两个误差限ε1和ε2,这样的参数是可以人为设定的。可以由opts=optimset命令得出方程求解的控制选项opts,该变量是一个结构体型的数据,其中很多成员变量是可以修改的,例如,AbsTol是绝对误差限,与前面介绍的常数ε1是一致的,更常用的是相对误差限RelTol,另外,函数值误差限FunTol与常数ε2是一致的,所以可以通过设置这些误差限来控制求解的精度。
除了这两个选项之外,迭代步数MaxIter和方程函数调用次数MaxFunEvals也经常需要重新设置,否则,在求解方程过程完成之后可能会给出提示,指明求解次数超限。在这种情况下,一方面可以将这两个选项设置为更大的值,另一方面,也可以将得出的结果作为初值,重新调用fsolve()函数继续求解,直至找出方程的解为止。这样的方法还可以与循环结构配合使用。
下面将通过例子演示求解精度的设定与控制方法。
例2-18 试选择幺矩阵作初值,重新求解例2-17中的Riccati代数方程。
解 如果重新设置求解精度控制选项ff,在调用语句中可以直接使用该控制变量,得出更精确的解,这时,误差矩阵的范数为2.5020×10−15,比默认精度有了明显的提高。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P34_29689.jpg?sign=1738892835-p2wU5PryUgi1dyJzKmXJMo3QvLZx1oBZ-0-3ef2661179baa076e517638e0c03e4f5)
例2-19 试求解下面的Riccati代数方程。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P34_29691.jpg?sign=1738892835-lA1NiwZfN9avc5FZcoLOt2ofqNbi0EWS-0-71a0620c3bd2c76a7e0b4f1650b16a3f)
解 由下面的语句可以尝试求解Riccati方程。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P35_29692.jpg?sign=1738892835-BHGhByEYhGrKwfCOOB76Od8AfdNj4izo-0-47ccc14388aaad464114d34fe37fa782)
不过求解过程中会得出“No solution:(A,B)may be uncontrollable or no solution exists”((A,B)可能不可控或方程无解)错误信息。尽管are()函数失效,仍然可以尝试求解原方程的数值解方法。例如,可以由下面的语句直接求解。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P35_29694.jpg?sign=1738892835-iZKMXX4qHvksCzIvexm6WmVeo8GF0WYx-0-647ba40030dcadea34d04b044e143f8d)
得出的一个解如下,该解的误差范数为1.2515×10−14。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P35_29697.jpg?sign=1738892835-n8u9BoLXZNuMBWVLBFoZNtfmBWiqFcun-0-a8a38663937e64667e1a8d8ea1c7cec2)
2.3.4 方程的复域求解
使用fsolve()的另一个优势是,当初值选作复数时,有可能得出方程的复数根。另外,该函数还可以求取复数系数方程的数值解。
例2-20 试求例2-17方程的复数根。
解 如果选择复数初值,则也有可能得出方程的复数根,同时可以验证,该根的共轭复数矩阵也满足原始方程。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P35_29699.jpg?sign=1738892835-aaBqopX4ySe2FXUcrITfKvWmJ3rtNgZy-0-ca9c33cd4d19f3d31e9e0596dd6990fd)
由该初值可以搜索出的解如下,相应的误差为1.1928×10−14。对这个例子还可以同步求出方程根的共轭矩阵,经检验该矩阵也满足于原方程。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P35_29701.jpg?sign=1738892835-YCPNbsz2llTbYbFjnQaFKWtKNoBKhdaw-0-9d4ffbba3e8008b7ae1aeb42b944f8a9)
对这个具体问题而言,如果初始值选择成实部为幺矩阵,虚部为单位阵的矩阵,则搜索出的将是实数解,这也说明由复数初始搜索点出发也能搜索出实数根。不过值得注意的是,这样搜索出的结果可能带有微小的虚部,其范数在该例子中为6.4366×10−19,应该由real()函数提取出来。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P35_29703.jpg?sign=1738892835-ttTKJDYab95D89448OEqyGRaBEy580ZF-0-48c69bec9e1c6fd399e0f3a735e0dd28)
例2-21 如果Riccati代数方程的系数矩阵A变成复数矩阵,试重新求解该方程。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P35_29705.jpg?sign=1738892835-1d3cFkg0VMvVo5MVpS3MstV5oagNwe3r-0-686b4d492c542c832075a391b276a761)
解 可以将各个矩阵输入到MATLAB的工作空间,然后使用匿名函数描述原方程。注意,原方程使用的是矩阵A的直接转置AT,不是Hermite转置AH,所以在匿名函数中应该使用A.',而不能使用A',否则方程描述是错误的,求解就没有意义了。可以使用下面命令直接求解并检验原方程,并求出解的共轭复数矩阵。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P36_29707.jpg?sign=1738892835-JxQ1fkKMVux38ICDtYhpajqc0ovTtMo8-0-26daa155794104e0f5a984b1ada68633)
这样可以得出方程的解为
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P36_29708.jpg?sign=1738892835-xa3OD0JVecL785Cvp1TMmDJhYnDw5lML-0-6e5f18610def5d20aa4a4d1920d1c060)
如果将解代入原方程,则得出误差矩阵的范数为5.4565×10−13。对这个特定的例子而言,即使使用实数起始搜索矩阵,也能搜索出方程的复数解。此外,虽然能直接求出解矩阵的共轭复数矩阵,但显然该矩阵不满足原方程,这与实数矩阵方程不同,因为在实数矩阵方程中,如果找到了方程的一个解,其共轭矩阵一般会满足原方程。
如果采用MATLAB控制系统工具箱的are()函数求解方程的解:
>> X=are(A,B,C), norm(f(X))
得出矩阵方程的解如下,在求解过程中也没有任何警告或错误信息,不过试图将该解代入原方程后得出的误差过大,为10.5210,说明得出的解根本不满足原始方程。
![](https://epubservercos.yuewen.com/2E5615/16499866905000206/epubprivate/OEBPS/Images/Figure-P36_29711.jpg?sign=1738892835-EEk4utsVpFExkmJbOwCNP1x8WvY9wsZO-0-aa0cd0da03e1eab8253e4b35cd5ca23b)