2.5 多解矩阵方程的求解

尽管前面介绍的vpasolve()函数能够一次性求出某些方程全部的解，但对一般矩阵方程，尤其是非线性矩阵方程而言是无能为力的，也没有其他MATLAB程序能够求解任意的非线性矩阵方程。

前面介绍了由给定初值求解函数的几种方法，但这些方法很难一次性求解多解非线性方程所有的解，所以应该构建更简单的函数完成这样的任务，本节给出一种求解的思路，并依照该思路编写出MATLAB通用程序，试图得出方程感兴趣区域内全部的解，并再扩展一步该方法，试图得出方程的全部准解析解。

2.5.1 方程求解思路与一般求解函数

由前面给出的求解函数可见，如果选定了一个初值，则可以通过随机数矩阵生成的方式产生一个初始搜索点，得出方程的解。更一般地，可以建立一个循环结构实现这样的操作。由初始搜索点，调用fsolve()函数得出方程的一个解，如果这个解已经被记录，则可以比较这个解和已记录解的精度，如果新的解更精确，则用这个解取代已记录的解，否则舍弃。如果这个解是新的解，则记录该解。

如果这个循环结构设计成死循环，则有望得出方程全部的解。根据这样的思路可以编写出一个通用的求解函数，其实这个函数以前曾经发布了多个版本，这个版本中，特别增加了一些处理，例如，可以尝试零矩阵是不是方程的孤立解；如果找到的根比以前存储的更精确，则替换该根；如果找到的新解为复数，则检验其共轭复数是不是方程的根。基于这样的考虑，编写了下面的求解函数。这个函数有一个特点，运行的时间越长，得出的结果可能越精确。

该函数调用底层公用函数default_vals()，其内容在其他卷中有过描述，为方便起见，这里重新给出。

方程求解函数more_sols()的调用格式为

    more_sols(f,X0,a,ϵ,tlim,opts)

式中，f为方程的函数句柄，可以由匿名函数与M函数描述原代数方程；X0为三维数组，用于描述解的初值，如果首次求解方程，建议将其设置为zeros(n,m,0)，即空白三维数组，n和m为解矩阵的维数；方程的解被自动存在MATLAB工作空间中的三维数组X中，如果想继续搜索方程的解，则应该在X0的位置填写X；a的默认值为1000，表示在[−500,500]区间内大范围搜索方程的解；ϵ的默认值为eps；tlim的默认值为30，表示30s没有找到新的解就自动终止程序；还可以指定求解的控制选项opts，默认值为optimset。a还可以取为复数，表示需要求取方程的复数根。另外，a还可以给定为求解区间[a,b]。

例2-36 试重新求解例2-10中的一元超越方程。

解 从图2-5给出的曲线可见，该方程在期望的区域内有6个交点，利用编写的more_sols()函数可以求出这6个交点，并在图2-11上直接标注出来。得出方程的解为x1=[1.4720,4.6349,7.7990,10.9601,14.1159,17.2666]。

例2-37 试求解例2-11给出的代数方程在−2π≤x1,x2≤2π区域内全部的根。

解 利用下面自然的方式可以直接求解非线性代数方程，先用匿名函数的形式描述联立方程，这样就可以调用more_sols()函数求取方程在感兴趣区域内全部的数值解。这里将A可以选作13，比4π略大的值，得出的解个数多于在感兴趣区域内的解。

可以提取出感兴趣区域内全部的根，可以发现，总的根的个数为110个，可以将得到的解叠印到图解法的图形上，如图2-12所示。

例2-38 试用数值方法重新求解例2-32中的广义Riccati方程。

解 和以前一样，先输入几个矩阵，然后由匿名函数的方式描述方程，再给出求解命令就可以直接求解方程，得到方程的8个实数根。

继续求解方程则可以得出方程全部的复数根，如果将A设置成复数量。这样总共可以得到方程的20个根，与例2-32得出的结果完全一致。

    >> more_sols(F,X,1000+1000i)

例2-39 如果例2-32中的广义Riccati方程变换成DX+XA−BX2+C=0，试用不同的方法求解该方程。

解 如果用more_sols()函数直接求解，则可以得到19个实数根。

如果在复数范围内搜索方程的根，总共可以找到57个根。现在尝试vpasolve()函数的准解析解求解方法，经过5873.8s的等待，可以得到60个根。

两种方法得出的根的数目差不多，不过仔细对比得出的根，如x11的值，可以发现，由后者得出的根有11个位于无穷远处，其幅值大于1010，而前者得出的都是在合理范围内的数值。两种方法得出的x11对比在表2-1中给出。此外，经检验，x11=38的根是第21个根，将其提取并代入原方程，得出误差矩阵的范数为18981334.363，不满足原方程。x11=−1.5146的根是第15个根，也不满足原始方程。由此可见，采用vpasolve()函数求解，即使处理多项式矩阵方程，得出的结果也是不可靠的。

例2-40 试求解非线性矩阵方程。

eAXsinBX−CX+D=0

其中A、B、C和D矩阵在例2-32中给出，试求出该方程的全部实根。

解 可以用下面的语句直接求解这里给出的复杂非线性矩阵方程，已经找到122个实根。用户还可以自己尝试，看看能不能得出更多的实根（根的存储文件为data2_36.mat）。

2.5.2 伪多项式方程的求解

伪多项式方程是非线性矩阵方程的一个特例，因为未知矩阵实际上是一个标量。这里将给出伪多项式方程的定义，并给出求解方法。

定义2-10 伪多项式（pseudo-polynomial）方程的一般数学形式为

式中，αi为实数。

可见，伪多项式方程是常规多项式方程的扩展，求解方法可能远难于普通多项式方程的求解方法。这里将探讨不同方法的可行性。

例2-41 试求解伪多项式方程[5]x2.3+5x1.6+6x1.3−5x0.4+7=0。

解 一种容易想到的方法是引入新变量z=x0.1，这样原方程可以映射成关于z的多项式方程如下：

z23+5z16+6z13−5z4+7=0

可以求出，该方程有23个根，再用x=z10就可以求出方程全部的根。这样的思想可以由下面的MATLAB语句实现。

不过，把这样得出的解代回到原来的方程可以发现，绝大部分的根都不满足原有的伪多项式方程。原方程到底有多少个根呢？上面得到的x只有两个根满足原方程，即x=−0.1076±0.5562j，其余的21个根都是增根。由下面语句也可以得出同样两个根。

从数学角度看，这对真实的根位于第一Riemann叶上。其余的解位于其他Riemann叶上，都是方程的增根，但不满足原方程。

例2-42 试求解无理阶次伪多项式方程。

解 由于这里的阶次是无理数，无法将其转换成前面介绍的普通多项式方程，所以more_sols()函数就成了求解这类方程的唯一方法，可以由下面的语句直接求解该方程。可见，该无理阶伪多项式方程只有两个根，位置为s=−0.0812±0.2880j。

将得出的解代回原方程，则可见误差为9.1551×10−16，该解在数值意义下足够精确。从这个例子还可以看出，即使阶次变成了无理数，并未给求解过程增加任何麻烦，求解过程和计算复杂度与前面的例子完全一致。

2.5.3 高精度求解函数

在这里给出的more_sols()函数中，核心求解工具是fsolve()函数，若将其替换为高精度的vpasolve()函数，则可以编写出高精度的非线性函数求解程序。

该函数的调用格式为more_vpasols(f,X0,A,tlim)，其中还嵌入了为其设计的底层支持子函数sol2vec()，将得出的解转换成行向量。输入变量f可以为符号型的行向量来描述联立方程，初始矩阵X0指定为zeros(0,n)，其中，n为未知数的个数。其他的输入变元与前面介绍的more_sols()函数是一致的。返回的变量X(i,:)存储找到的第i个解。值得指出的是，more_vdpsols()函数的运行速度比more_sols()函数慢得多，但精度也高得多。

例2-43 考虑例2-11中的联立方程，试找出−2π≤x,y≤2π范围内所有准解析解。

解 可以用下面的命令直接求解联立方程：

要检验得出方程根的精度，则首先应该提取出感兴趣区域内的根x0和y0，并对其进行排序，得出的根代入原方程后的误差范数为7.79×10−32，比例2-37中得出的精度要高得多，所需的时间在半个小时左右，也远高于more_sols()函数，用这样的方法只找到区域内的105个根。得出的图形类似于图2-12，不过有几个点缺失。

例2-44 试求解例2-40的高精度数值解。

解 可以试用下面语句求解方程，不过在用符号表达式描述方程时，并不能计算出方程的表达式f，所以不能调用后续的more_vpasols()函数，不能得出方程的高精度数值解，求解这类方程只能采用例2-40的普通数值解方法。