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模式。
满秩方程
假设有如下满秩方程:
用矩阵表示:
从几何意义来说,这种方程组的解是空间中三个平面的交点,为了更加直观展现,下面用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)
这个语句来 声明变量名称和个数。
欠秩方程
假设欠秩方程组如下:
此时仍可以用矩阵求逆的方法来做,只不过现在是求伪逆。按照线性代数的方法,我们知道欠秩方程组的解包含通解+特解,在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.])
超定方程
对于如下超定方程:
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,具体的可以查看官方帮助
总结了三种线性方程组的求解方法,以后想到了新的会增添内容,内容就先这么多了。