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的位置为
式中，f′(x0)为f(x)关于x的导函数在x0点的值。再由x1出发作切线，则可以得出x2，以x2出发找到x3，…。若已经找到了xk，则从该点出发可以搜索到下一个点。
如果|xk+1−xk|≤ε1或|f(xk+1)|≤ε2，其中，ε1与ε2为预先选定的误差容限，则可以认为xk为原方程的一个解。
图2-8 Newton–Raphson迭代法示意图
定义2-7 对多元向量函数f(x)，其Jacobi矩阵的定义为
定理2-3 更一般地，对多元方程f(x)=0，其中x为向量或矩阵，而非线性函数f(x)也是同维数的函数，则可以由下式迭代求出：
如果||xk+1−xk||≤ε1或||f(xk+1)||<ε2，则xk为方程的根。在定义中使用了Jacobi矩阵。
根据给出的算法，编写出MATLAB的通用求解函数：
该函数需要用户提供方程函数句柄f及其Jacobi矩阵的句柄d，还需要给出初值的列向量x0，并给出误差限ε。如果给出key，并令其值为1，则可以将中间搜索结果由x矩阵返回。
对单变量函数而言，如果需要提供给定函数的导函数本身就是个比较麻烦的事，使用可以考虑采用正割方法来近似函数的导数，搜索方程的根。这时，求解式（2-3-2）可以改写成
如果采用点运算，则这里的正割函数同样适用于多元方程的求解。
例2-13 选择初值x0=10，并求出例2-10中一元超越方程的一个根。
解 先由符号运算的方式推导出给定函数的一阶导数。
>> syms x; f=exp(-0.2*x)*sin(3*x+2)+cos(x), diff(f)
可以得出函数的导数为
有了函数及其导数，则可以用匿名函数表示它们，再设定搜索的初值为x0=10，这时调用Newton–Raphson求解算法，则可以得出方程的解搜索的中间点为[10,10.8809,11.0700,11.0593,11.0593]T，方程的解为x=11.0593，将其代入原方程，则可以得出误差为−5.1348×10−16。可以看出，利用这样的方法求解精度还是比较高的。整个求解过程的示意图如图2-9所示，可以看出经过几步迭代就可以得出方程的解。
图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算法，但其优势是无须用户提供函数的导函数，所以该算法是有意义的。
图2-10 一元方程求解的中间过程
例2-15 试求解例2-11中的二元方程，以(1,1)点为初值搜索到一个解。
解 由于同时含有自变量x和y，这样的方程是不能直接求解的，需要将其改写成x的方程，最简单的方法是令x1=x，x2=y，这样，原始的方程可以改写成
Jacobi矩阵不易用手工的方法推导出来，是需要通过解析推导求解的，所以应该在符号运算框架下输入原始函数，并通过jacobian()函数计算出该矩阵。
可以推导出函数的Jacobi矩阵为
有了原函数与Jacobi矩阵，则可以手工写出这两个函数的匿名函数，然后调用基于Newton–Raphson算法的求解函数。
由该初值点出发得出的方程的解为x=[5.1236,−12.2632]T，中间点的个数为18。代入原方程后得出的误差矩阵范数为3.9323×10−12。从求解的结果看，确实通过该函数可以求出方程的一个根，不过从实际操作看，这样的方法未免过于复杂，由于需要推导Jacobi矩阵，并将其矩阵用匿名函数的形式手工表示出来，该过程比较容易出错。
若采用正割方法，则可以给出下面的语句，得出方程的根为x=[1.0825,−1.1737]T，代入方程得出误差为2.2841×10−14。
虽然这时函数调用无须用户提供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中的结果相仿。不过该函数的优势是无须用户提供方程的导函数，使用该函数更适合于实际应用，建议采用该方法直接求解方程。
如果将初值修改为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代数方程的数学形式为
式中，每个矩阵都是n×n矩阵。
Riccati方程是以意大利数学家Jacopo Francesco Riccati（1676−1754）命名的，本意是对应于一类一阶微分方程，其原型要求B为正定矩阵，C为对称矩阵。后来因为微分方程难以求解，所以将其简化成上面的Riccati代数方程。
从数学角度看，这几个矩阵可以为任意矩阵。在控制科学领域，可以考虑采用控制系统工具箱提供的are()函数直接求取方程的数值解，不过该函数只能求出Riccati方程的一个根。如果想获得方程全部的根，则可以考虑调用vpasolve()函数直接求解。下面将给出例子演示Riccati方程与多项式矩阵方程的方法。
例2-17 试求解下面的Riccati代数方程。
解 可以先输入这几个矩阵，然后调用控制系统工具箱中的are()函数求解Riccati方程，并计算得出的误差矩阵范数。
由控制系统工具箱的are()函数可以直接得出方程的一个根如下，代入原方程可见，该解导致的误差为1.3980×10−14，说明该解的精度是比较高的。
由于该方程是二次型方程，人们很自然就可以想到，该方程可能有多个根，如何求出其他的根呢？显然可以尝试选择一个不同的初始矩阵，例如幺矩阵，从该矩阵求解方程的根。
得出方程的解如下，对应的误差为6.3513×10−9。
显然，这时得出的解与are()函数得出的解是不一致的。该函数返回的f1矩阵就是方程解处的函数值，即f(X1)。另外，在这个例子中返回的flag值为1，因为它为正数，所以表示求解成功。另外返回的cc信息包括如下内容，表明迭代步数为11，函数调用次数为102，说明求解效率还是很高的。
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，比默认精度有了明显的提高。
例2-19 试求解下面的Riccati代数方程。
解 由下面的语句可以尝试求解Riccati方程。
不过求解过程中会得出“No solution:(A,B)may be uncontrollable or no solution exists”（(A,B)可能不可控或方程无解）错误信息。尽管are()函数失效，仍然可以尝试求解原方程的数值解方法。例如，可以由下面的语句直接求解。
得出的一个解如下，该解的误差范数为1.2515×10−14。
2.3.4 方程的复域求解
使用fsolve()的另一个优势是，当初值选作复数时，有可能得出方程的复数根。另外，该函数还可以求取复数系数方程的数值解。
例2-20 试求例2-17方程的复数根。
解 如果选择复数初值，则也有可能得出方程的复数根，同时可以验证，该根的共轭复数矩阵也满足原始方程。
由该初值可以搜索出的解如下，相应的误差为1.1928×10−14。对这个例子还可以同步求出方程根的共轭矩阵，经检验该矩阵也满足于原方程。
对这个具体问题而言，如果初始值选择成实部为幺矩阵，虚部为单位阵的矩阵，则搜索出的将是实数解，这也说明由复数初始搜索点出发也能搜索出实数根。不过值得注意的是，这样搜索出的结果可能带有微小的虚部，其范数在该例子中为6.4366×10−19，应该由real()函数提取出来。
例2-21 如果Riccati代数方程的系数矩阵A变成复数矩阵，试重新求解该方程。
解 可以将各个矩阵输入到MATLAB的工作空间，然后使用匿名函数描述原方程。注意，原方程使用的是矩阵A的直接转置AT，不是Hermite转置AH，所以在匿名函数中应该使用A.'，而不能使用A'，否则方程描述是错误的，求解就没有意义了。可以使用下面命令直接求解并检验原方程，并求出解的共轭复数矩阵。
这样可以得出方程的解为
如果将解代入原方程，则得出误差矩阵的范数为5.4565×10−13。对这个特定的例子而言，即使使用实数起始搜索矩阵，也能搜索出方程的复数解。此外，虽然能直接求出解矩阵的共轭复数矩阵，但显然该矩阵不满足原方程，这与实数矩阵方程不同，因为在实数矩阵方程中，如果找到了方程的一个解，其共轭矩阵一般会满足原方程。
如果采用MATLAB控制系统工具箱的are()函数求解方程的解：
>> X=are(A,B,C), norm(f(X))
得出矩阵方程的解如下，在求解过程中也没有任何警告或错误信息，不过试图将该解代入原方程后得出的误差过大，为10.5210，说明得出的解根本不满足原始方程。
