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

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

满秩方程

假设有如下满秩方程：

⎧⎩⎨3x+4y−2z=74x+y+3z=6x+y+7z=5

用矩阵表示：

⎡⎣⎢341411−237⎤⎦⎥⎡⎣⎢xyz⎤⎦⎥=⎡⎣⎢765⎤⎦⎥

从几何意义来说，这种方程组的解是空间中三个平面的交点，为了更加直观展现，下面用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 onmesh(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 npimport mayavi.mlab as mlimport matplotlib.pyplot as pltx,y=np.mgrid[-5:5:0.2,-4:4:0.1] #生成网格z1 = (3*x+4*y-7)/2z2 = (6-4*x-y)/3z3 = (5-x-y)/7ml.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 splotimport sympy as spx,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) %符号解

求解结果如下：

>> equationsx1 =    0.8723    1.2979    0.4043x2 =    0.8723    1.2979    0.4043x3 =    0.8723    1.2979    0.4043x4 = 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 npa = 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]: xOut[82]: array([[ 0.87234043],       [ 1.29787234],       [ 0.40425532]])

Sympy:

import sympy as spa = sp.Matrix([[3,4,-2],[4,1,3],[1,1,7]])b = sp.Matrix([[7],[6],[5]])x = sp.symarray(x,3)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) 这个语句来 声明变量名称和个数。

欠秩方程

假设欠秩方程组如下：

{3x+4y−2z=74x+y+3z=6

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

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

RESULTS:

x1 = -14/13  17/13      1Warning: System is rank deficient. Solution is not unique. x2 = 17/13 10/13     0x = 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 npa = np.array([[3,4,-2],[4,1,3]])b = np.array([[7],[6]])x = np.linalg.pinv(a).dot(b)

RESULTS:

In [138]: xOut[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.])

超定方程

对于如下超定方程：

⎧⎩⎨⎪⎪⎪⎪3x+4y−2z=74x+y+3z=6x+y+7z=52x−y+3z=4

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\bx2 = pinv(a)*b

RESULTS:

x1 =    1.2367    0.8868    0.3999x2 =    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]: x1Out[144]: (array([[ 1.23667528],       [ 0.88682789],       [ 0.39985912]]), array([ 2.8410425]), 3, array([ 8.87798703,  5.91668195,  2.48479798]))In [145]: x2Out[145]: array([[ 1.23667528],       [ 0.88682789],       [ 0.39985912]])

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

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