线性方程组求解——基于MTALAB/Octave,Numpy,Sympy和Maxima

---

2018/06/14更正
sympy代码运行出错，满秩的情况只要修改代码

x = sp.symarray(x,3)

为

x = sp.symarray('x',(3,1))

线性代数里一个重要的内容就是线性方程的求解，解方程其实从我们初中的时候就已经接触了，这篇文章记录的是对满秩方程（恰定方程）、欠秩方程（欠定方程）和超定方程三种线性方程的计算机求解方法，使用了MATLAB/Octave,Numpy,Sympy和Maxima来实现(有些可能是只是其中的几种)，除了MATLAB外，其他都是开源免费的，Octave和matlab最相似，大部分语法兼容，Numpy和Sympy是基于IPython平台，我自己是用轻量级Python IDE——IEP3.7(Pyzo老版本)，该IDE支持IPython模式。

满秩方程

假设有如下满秩方程：

$$\begin{cases} 3x+4y-2z=7 \\ 4x+y+3z=6\\ x+y+7z=5 \end{cases}$$
用矩阵表示：

$$\left[ \begin{array}{ccc} 3&4&-2\\ 4&1&3\\ 1&1&7 \end{array} \right] \left[ \begin{array}{ccc} x\\ y\\ z \end{array} \right] = \left[ \begin{array}{ccc} 7\\ 6\\ 5 \end{array} \right]$$
从几何意义来说，这种方程组的解是空间中三个平面的交点，为了更加直观展现，下面用Maxima,Mayavi和Sympy绘制上述方程组每个方程代表的平面。

MATLAB/Octave

[x,y]=meshgrid (-5:0.1:5);
z1=(3.*x+4.*y-7)/2;
z2=(6-4.*x-y)/3;
z3=(5-x-y)/7;

mesh(x,y,z1)
hold on
mesh(x,y,z2)
mesh(x,y,z3)
hold off

MATLAB绘制出的图形：

Maxiam:

wxplot3d([(3*x+4*y-7)/2,(6-4*x-y)/3,(5-x-y)/7, [x,-5,5], [y,-5,5]]);
或者：
plot3d([(3*x+4*y-7)/2,(6-4*x-y)/3,(5-x-y)/7, [x,-5,5], [y,-5,5]],[plot_format,xmaxima]);
或者：
plot3d([(3*x+4*y-7)/2,(6-4*x-y)/3,(5-x-y)/7, [x,-5,5], [y,-5,5]],[plot_format,gnuplot])$

下图是gnuplot格式绘制的

下图是xmaxima格式绘制的
Mayavi:

import numpy as np
import mayavi.mlab as ml
import matplotlib.pyplot as plt

x,y=np.mgrid[-5:5:0.2,-4:4:0.1] #生成网格

z1 = (3*x+4*y-7)/2
z2 = (6-4*x-y)/3
z3 = (5-x-y)/7

ml.surf(x,y,z1,colormap="magma")    #colormap参数值可以自己设定
ml.surf(x,y,z2,colormap="spring")
ml.surf(x,y,z3,colormap="winter")
ml.show()

Sympy：

import sympy.plotting as splot
import sympy as sp

x,y,z=sp.symbols('x,y,z')
splot.plot3d((3*x+4*y-7)/2,(6-4*x-y)/3,(5-x-y)/7,(x,-5,5),(y,-5,5))

有了图形的直观理解后，下面来求解方程，首先用matlab来实现，是最简单的。
MATLAB/Octave：

a = [3,4,-2;4,1,3;1,1,7];
b = [7;6;5];
x1 = a\b        % 形如A*X=B形式的用左除
x2 = inv(a)*b   % 先求逆
x3 = pinv(a)*b  % 求伪逆
x4 = sym(a)\sym(b) %符号解

求解结果如下：

>> equations
x1 =
    0.8723
    1.2979
    0.4043
x2 =
    0.8723
    1.2979
    0.4043
x3 =
    0.8723
    1.2979
    0.4043
x4 =
 41/47
 61/47
 19/47

Maxima:

solve([3*x+4*y-2*z-7,4*x+y+3*z-6,x+y+7*z-5],[x,y,z]);

Results:

[[x=41/47,y=61/47,z=19/47]]

Numpy:

import numpy as np
a = np.array([[3,4,-2],[4,1,3],[1,1,7]])
b = np.array([[7],[6],[5]])
x = np.linalg.solve(a,b)

Results:

In [82]: x
Out[82]: 
array([[ 0.87234043],
       [ 1.29787234],
       [ 0.40425532]])

Sympy:

import sympy as sp
a = sp.Matrix([[3,4,-2],[4,1,3],[1,1,7]])
b = sp.Matrix([[7],[6],[5]])

x = sp.symarray('x',(3,1))
sp.solve(a*x-b)

Results:

In [87]: sp.solve(a*x-b)
Out[87]: 
{[[ 0.87234043]
  [ 1.29787234]
  [ 0.40425532]]_0: 41/47, [[ 0.87234043]
  [ 1.29787234]
  [ 0.40425532]]_1: 61/47, [[ 0.87234043]
  [ 1.29787234]
  [ 0.40425532]]_2: 19/47}

由结果来看可以知道4种工具求解的函数语法基本一致，matlab本来就是为矩阵运算打造的语言，计算最方便，maxima适合方程组数和变量较少的情况，当然，博主 本人还没怎么用maxima处理矩阵，可能也有函数求解矩阵方程，暂时不清楚；sympy给出的是符号解和数值解，对于求解的变量需要增加x=symarray(x,3) 这个语句来 声明变量名称和个数。

欠秩方程

假设欠秩方程组如下：

$$\begin{cases} 3x+4y-2z=7 \\ 4x+y+3z=6 \end{cases}$$

此时仍可以用矩阵求逆的方法来做，只不过现在是求伪逆。按照线性代数的方法，我们知道欠秩方程组的解包含通解+特解,在matlab中可以用null()函数求得通解，用左除求得特解。
MATLAB/Octave:

a = [3,4,-2;4,1,3];
b = [7;6];
syms k
x1 = null(sym(a))
x2 = sym(a)\sym(b)
x = k*x1+x2

RESULTS:

x1 =
 -14/13
  17/13
      1
Warning: System is rank deficient. Solution is not unique. 
x2 =
 17/13
 10/13
     0
x =
 17/13 - (14*k)/13
 (17*k)/13 + 10/13
                 k

Maxima:

solve([3*x+4*y-2*z-7,4*x+y+3*z-6],[x,y,z]);

RESULTS:

[[x=-(14*%r2-17)/13,y=(17*%r2+10)/13,z=%r2]]

Sympy:

a=sp.Matrix([[3,4,-2],[4,1,3]])
b=sp.Matrix([[7],[6]])
x,y,z=sp.symbols('x,y,z')
x=sp.symarray(x,3)
t = sp.solve(a*x-b)

RESULTS:

Out[134]: {x_1: 17*x_2/13 + 10/13, x_0: -14*x_2/13 + 17/13}
In [135]: a.nullspace()
Out[135]: 
[Matrix([
 [-14/13],
 [ 17/13],
 [     1]])]

Numpy:
Numpy没有nullspace()这样的函数，我网上查了一下可以用已有的函数自己实现，但因为自己当时学线性代数对SVD分解等不了解，就没有自己写代码了，有兴趣的同学可以自己写一下，在numpy解欠秩方程用的是伪逆pinv()，其实这个在本科线性代数也是没有学的，无意中了解到的，但是很好用，因为伪逆也是由SVD分解之后转换而来，伪逆适用于所有类型的矩阵，不像一般的逆矩阵只适用于方形矩阵（square matrix）

import numpy as np
a = np.array([[3,4,-2],[4,1,3]])
b = np.array([[7],[6]])
x = np.linalg.pinv(a).dot(b)

RESULTS:

In [138]: x
Out[138]: 
array([[ 1.19571865],
       [ 0.90519878],
       [ 0.10397554]])

验证：

In [141]: a[0,:].dot(x)
Out[141]: array([ 7.])

In [142]: a[1,:].dot(x)
Out[142]: array([ 6.])

超定方程

对于如下超定方程：

$$\begin{cases} 3x+4y-2z=7 \\ 4x+y+3z=6 \\ x+y+7z=5\\ 2x-y+3z=4 \end{cases}$$

matlab用左除”\”返回基于最小二乘法计算出的解，maxima还不知，numpy提供了lstsq()函数和伪逆的方法求解，sympy暂时不知。

MATLAB/Octave:

a=[3,4,-2;4,1,3;1,1,7;2,-1,3];
b=[7;6;5;4];
x1 = a\b
x2 = pinv(a)*b

RESULTS:

x1 =
    1.2367
    0.8868
    0.3999
x2 =
    1.2367
    0.8868
    0.3999

Numpy:

a = np.array([[3,4,-2],[4,1,3],[1,1,7],[2,-1,3]])
b = np.array([[7],[6],[5],[4]])

x1 = np.linalg.lstsq(a,b)
x2 = np.linalg.pinv(a).dot(b)

RESULTS:

In [144]: x1
Out[144]: 
(array([[ 1.23667528],
       [ 0.88682789],
       [ 0.39985912]]), array([ 2.8410425]), 3, array([ 8.87798703,  5.91668195,  2.48479798]))

In [145]: x2
Out[145]: 
array([[ 1.23667528],
       [ 0.88682789],
       [ 0.39985912]])

返回的是tuple类型，第一个元素是基于最小二乘法求出的解，第二个数residues，第三个是rank，具体的可以查看官方帮助

总结了三种线性方程组的求解方法，以后想到了新的会增添内容，内容就先这么多了。