2.4 联立方程组的精确求解

前面介绍过，利用图解方法只能求出给定方程的实数根，并不能求出方程的复数根，具体例子可以参见例2-12。另外，如果联立方程有多个实数根，则只能用图形方法绘制出根所在的位置，并不能直接得出根的具体值，需要逐个根进行局部放大求解，求解过程比较烦琐。此外，前面介绍的数值求解方法由于每次只能求出方程的一个根，使用起来有时也不方便。

MATLAB工具箱提供了代数方程的解析求解函数solve()，可以直接用于求解解析解存在的代数方程。如果方程的解析解不存在，则可以采用vpasolve()函数求取方程的高精度数值解，解的误差可能达到10−30或更高精度，远大于双精度数据结构下的数值解。本书称这类解为准解析解，以区别于方程的解析解与双精度意义下的数值解。

2.4.1 低阶多项式方程的解析求解

一元一次和一元二次方程可以利用solve()函数直接求解，该函数还可以用于含有其他参数的方程求解。不过如果有其他参数的存在，则可能利用现有的solve()函数不易得出三次或四次方程的解析解，尽管这些方程是存在代数解法的，除非给出特别的控制选项，后面将通过例子演示。

solve()函数的调用格式为

式中，待求解的方程由符号表达式eqni表示，自变量由xi表示。返回的解是一个结构体型变量，其解由S.xi直接提取。在调用格式中eqni可以是单个的方程也可以是向量、矩阵描述的一组方程，还可以将所有的方程描述成一个向量与矩阵符号表达式eqn1，直接求解这些方程。

当然，用下面的调用格式还可以直接获得方程的解。

[x1,x2,…,xn]=solve(...),%输入变量表示与前面一致

从函数的调用和使用方面看，这种直接返回变量的调用格式比返回结构体变量的格式更实用，所以本书尽量使用这样的格式。

例2-22 试重新求解例2-2中的鸡兔同笼问题。

解 声明符号变量，并将方程用符号表达式表示出来，则可以调用solve()函数直接求解给定的方程。注意，方程中的等号应该由双等号表示。

得出方程的解为x1=23，x2=12。返回变量的名字还可以选择为其他的变量，此外，如果方程右边为零，则可以省略等号。例如，上面的求解语句还可以修改成

    >> syms x1 x2; [x0,y0]=solve(x1+x2-35,2*x1+4*x2-94)

例2-23 中国唐代数学家王孝通所著《缉古算经》中有一道应用题，翻译成现代数学符号可以写出下面的三次方程，试求解该方程。

解 利用符号运算方法可以直接求解这个三次方程。

得出方程的三个解如下。王孝通只得到了31这个解，由于它为整数，是原应用题的解，其他两个是负数，不是应用题的解。

例2-24 试求解下面的联立方程组，并给出方程有实数根的条件。

解 可以看出，如果对方程进行简单变换，则可以得到关于x或y一元二次方程，利用相应的求根公式，则可以写出方程的解析解。不过从给出的表达式看，如果想手工求解这样的方程并不是一件简单的事情，所以应该利用计算机来完成这样烦琐的推导。

先声明必要的符号变量，则可以将两个方程用符号表达式描述出来，再给出下面的直接求解命令。

则可以得出方程的一对解为

由得出的结果看，如果期望方程有实数根，则根号下的表达式应该大于等于零，即

c4−48ac−8a2c−12bc2−12a3−16c2−16ab+64≥0

例2-25 试在符号运算的框架下求取例2-7中复系数多项式方程的解析解。

解 由于已知方程的解析解存在，所以可以调用solve()函数解方程，最终得出方程的解析解为[−2,−1−j,−1−j,−1−j]T。

例2-26 试求解四次方程的解析解。

7s4+119s3+756s2+2128s+2240=0

解 由于四次方程有代数解求根公式，所以可以将其输入MATLAB环境，则可以调用solve()函数尝试求解。

得出方程的解为x=−5，−4，−4，−4。事实上，虽然四次方程有自己的代数解求根公式，得出的解也可能是无理数，换句话说，真正意义下的解析解也可能不存在。例如，如果原方程左侧加1，则方程仍然是四次多项式方程，这时方程的解析解可能不存在。

    >> f=f+1; x=solve(f)

这样的方程由solve()函数是不能直接求解的，因为直接得出的结果如下：

        root(z^4 + 17*z^3 + 108*z^2 + 304*z + 2241/7, z, 1)
        root(z^4 + 17*z^3 + 108*z^2 + 304*z + 2241/7, z, 2)
        root(z^4 + 17*z^3 + 108*z^2 + 304*z + 2241/7, z, 3)
        root(z^4 + 17*z^3 + 108*z^2 + 304*z + 2241/7, z, 4)

表明方程的解是z4+17z3+108z2+304z+2241/7=0的多项式方程的四个根。如果实在想得出方程的数值解，则可以调用vpa()函数或直接调用vpasolve()函数直接求解方程，得出的误差向量范数为4.7877×10−36。

    >> x1=vpa(x), norm(subs(f,s,x1))

其实，若真的想求出三次或四次方程的解析解，还是可以通过设定'MaxDegree'选项实现的，不过得出的解很冗长。下面将通过例子演示得出的解。

例2-27 试求出例2-26改变后四次方程的解析解。

解 可以重新求解该方程，由于是四次方程，所以将'MaxDegree'选项设置为4。

得出的结果还是很烦琐的，经过细心的手工变量替换与化简可得

式中，

有了这样的求解公式，则可以得出任意精度的方程的解。例如，由下面的语句可以得出误差为7.9004×10−104。

    >> norm(vpa(subs(f,s,x),100))

2.4.2 多项式型方程的准解析解

对于一般的多项式型代数方程而言，采用vpasolve()函数有望得出方程全部的准解析解，而一般非线性代数方程只能得到一个准解析解。本节将通过例子演示各类方程的准解析解方法。

例2-28 试求解例2-12中给出的联立方程。

解 由图解法并不能有效求解例子中的联立方程，因为原方程既含有实数根，也含有复数根。可以先将方程用符号表达式的方式输入给计算机，然后调用vpasolve()函数，则可以直接求解该方程组。

得出的解为

得出的结果在x0和y0向量中返回，所以要检验得出的误差并不是一件容易的事，因为要重新输入方程的表达式，再求值。另一种求解的方法是将两个方程用符号变量表示出来，再直接求解，这样得出的结果与前面是完全一致的。将解代入原方程，可以得出误差为1.0196×10−38。

从这里给出的例子看，求解这样复杂的高阶代数方程组的难度对用户而言，与求解鸡兔同笼问题是一样的，用户需要做的就是将方程用符号表达式表示出来送给求解函数，然后等待得出的解就行了。

更简单地，还可以用向量的形式表示两个方程，这样，再调用vpasolve()函数则可以重新求解相应的方程，得出的结果与前面语句完全一致。

例2-29 试求解下面看起来很复杂的代数方程。

解 这个方程与例2-12中给出的方程不一样，至少例2-12中的方程可以将一个方程代入另一个方程，最终手工得出一个高阶的多项式方程，而这个方程就没那么容易变换了。不用说求解方程，就是判定方程有多少个根，不借助计算机也不是一件容易的事。

不过不必考虑或担心这类底层问题，只需用下面的语句将原始的方程用规范的语句原原本本表示出来，就可以直接得出原方程的准解析解。

将得出的全部26个根代入原始方程，则能得出很小的计算误差，达到10−33级，说明该方程的各个解都是非常精确的。求解这样看起来难度极高的代数方程，对用户而言，求解的难度也与鸡兔同笼问题是一样的。

即使著名的Abel–Ruffini定理已经指明这类高阶多项式方程没有解析解，还是可以通过代数方程求解方法得出高精度的准解析解，得出解的精度是非常高的。

2.4.3 高次多项式矩阵方程的准解析解

定义2-8中描述了代数Riccati方程，该方程是关于X的二次型方程，是多项式型方程。如何用符号表达式表示Riccati方程是求解该方程的关键一步。这里先通过例子探讨使用vpasolve()函数求解该方程的方法。

例2-30 试求解例2-17中给出的Riccati代数方程全部的根。

解 如果想得出方程全部的根，则应该尝试vpasolve()函数。若想表示未知的X矩阵，则需要将其设置为符号变量，再构成符号矩阵，这样就可以由简单的符号表达式描述Riccati方程本身，再调用vpasolve()函数求解方程。

由符号运算可以推导出Riccati方程对应的联立方程为

经过23.04s的等待，可以得出方程全部的20组根，其中，第5、10、15、18这四组根为实数矩阵，其余为共轭复数矩阵。得出的四个实数矩阵分别为

例2-17中得出的实数根是原方程的第15个根，可以提取该根，代入方程检验，得出的误差为1.8441×10−39。

例2-31 试用准解析解方法求解例2-19中方程的全部根。

解 给出下面的命令即可求出方程的全部20个根，其中有8个实根，其余为共轭复数根，求解过程耗时为36.75s。

前面介绍了Riccati方程的求解方法，如果将Riccati方程中AT项替换成一个自由矩阵D，则可以引出广义Riccati方程。

定义2-9 广义Riccati代数方程的数学形式为

式中，各个矩阵都是n×n矩阵。

广义Riccati方程没有现有的MATLAB函数直接求解，仍可尝试vpasolve()函数直接求解，得出方程全部的根。

例2-32 若已知如下矩阵，试求解广义Riccati方程。

解 和以前一样，先输入几个矩阵，然后由符号表达式的方式描述方程，再给出求解命令就可以直接求解方程了。

经过25s的等待，由上述求解语句可以直接得出方程的20个根，其中，8个根为实数矩阵，其余的为共轭复数矩阵。

例2-33 本丛书卷III例5-42曾探讨了一个高次矩阵方程的求解问题。

AX3+X4D−X2BX+CX−I=0

其中，已知的矩阵为

试用这里给出的方法求取该方程的解矩阵。

解 很自然地，可以给出下面的求解命令，不幸的是，经过1843.8s的长时间等待，只得出了方程的一个根，并给出警告信息“Solutions might be lost”（根可能丢失了），说明用这样的方法也不能确保得出方程所有的根。

例2-34 试求解含有复数矩阵的Riccati方程的准解析解。

解 很自然地可以考虑使用下面的语句直接求解需要的复系数Riccati方程，同样的步骤和语句，由于涉及复系数，求解过程极其耗时，大约379s才能得出结果（实系数方程耗时23.04s，相差16倍）。可以得出方程的20个根，全部为复数根。

2.4.4 非线性代数方程的准解析解

如果用符号表达式可以描述出非线性联立方程组，则可以由vpasolve()函数直接求解方程，得出方程的解。与多项式类方程不同，这样得出的解很可能是众多解中的一个，如果原始方程含有多个解，则用户可以自己选择一个搜索的初值，从该初值x0出发直接搜索准解析解。相应的函数调用格式为

    x=vpasolve(eqn1,eqn2,···,eqnn,x1,x2,···,x,n,x0)

式中，x0为未知量向量的初始搜索点。如果原方程是多项式型的联立方程，则x0对求解结果没有影响。该函数的另一种调用格式为

    x=vpasolve(eqn1,eqn2,···,eqnn,x1,x2,···,xn,xm,xM)

式中，xm与xM为未知数向量x的下界与上界向量。

例2-35 试求解例2-11给出的代数方程的根，可以选择初始搜索点(2,2)。

解 利用下面自然的方式可以直接求解非线性代数方程。

得出方程的解如下，代入原方程则误差的范数为2.0755×10−37。由此可见，方程解的精度远高于前面介绍的数值解法精度。

x0=3.0029898235869693047992458712192

y0=−6.2769296697194948789764344182923

从给出的例子可见，尽管这里给出的方法能得出方程的高精度解，但仍然有一个问题尚未解决，就是如何一次性地求出如图2-6所示的所有交点坐标，即联立方程在感兴趣区域内所有的解，这将是2.5节需要解决的问题。