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模式。

满秩方程

假设有如下满秩方程：

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

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

⎡⎣⎢341411−237⎤⎦⎥⎡⎣⎢xyz⎤⎦⎥=⎡⎣⎢765⎤⎦⎥[34−2413117][xyz]=[765]

\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 pltx,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)/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 splot
import 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) %符号解

求解结果如下：

>> equations
x1 =0.87231.29790.4043
x2 =0.87231.29790.4043
x3 =0.87231.29790.4043
x4 =41/4761/4719/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) 这个语句来声明变量名称和个数。

欠秩方程

假设欠秩方程组如下：

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

\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/1317/131
Warning: System is rank deficient. Solution is not unique.
x2 =17/1310/130
x =17/13 - (14*k)/13(17*k)/13 + 10/13k

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.])

超定方程

对于如下超定方程：

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

\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.23670.88680.3999
x2 =1.23670.88680.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，具体的可以查看官方帮助

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

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

matlab中欠定方程组超定方程组_一篇文章入门大规模线性方程组求解
前面介绍过主要的线性方程组求解库,参考附录.求解大规模线性方程组是仿真软件求解器的底层技术,求解器时间基本都消耗在方程组求解上.线性方程组的解法比较成熟,方法也有很多,而且不同的方法对应不同类型方程组 ...
大型稀疏线性方程组求解技术——工业仿真的底层核心
背景在工业仿真领域,对各种现实世界的问题进行数值模拟时,如流体动力学分析.电磁场仿真.结构力学应力应变分析等,其控制方程通常是偏微分方程组,在经过不同方法的隐式离散之后最终都可转化为大型稀疏线性 ...
Pingouin: 基于pandas和numpy的统计包
Python网络爬虫与文本数据分析 pingouin是基于Pandas和numpy开发的Python3统计包.主要统计功能有方差分析多元线性回归中介效应分析卡方检验 Q-Q图贝叶斯因子信效 ...
22（线性方程组求解）高斯赛德尔迭代法
高斯赛德尔迭代法(线性方程组求解) [问题描述]为求解一个线性方程组,使用高斯赛德尔迭代法,采用欧几里得距离判定是否收敛.精度delta为1E-9,最大迭代次数为20. [输入形式]在屏幕上依次输入方 ...
大规模线性方程组求解
关于优化 1.大规模线性方程组求解 https://blog.tsingjyujing.com/game201/linear_eq_solver https://zhuanlan.zhihu.com/ ...
六、线性方程组求解--Jacobi和Gauss-Seidel迭代求解
六.线性方程组求解-Gauss消去.Jacobi和Gauss-Seidel迭代求解 1.消去法求解原理求解Ax = b,首先对增广矩阵(A|b)进行初等行变换为如下三种矩阵对角阵.上三角阵.下三角阵 ...
科学计算与数学建模-线性方程组求解的迭代法思维导图
第六章回归问题-线性方程组求解的迭代法
python中numpy数组的合并_基于Python中numpy数组的合并实例讲解
基于Python中numpy数组的合并实例讲解 Python中numpy数组的合并有很多方法,如 - np.append() - np.concatenate() - np.stack() - np. ...
python应用-scipy,numpy,sympy计算微积分
python应用-scipy,numpy,sympy计算微积分今天来讲一下使用python进行微积分运算,python有很多科学计算库都可以进行微积分运算,当然如果知晓微积分计算的原理也可以自己编程 ...
病态线性方程组求解（基于MATLAB）
A=[25 -300 1050 -1400 630; -300 4800 -18900 26880 -12600; 1050 -18900 79380 -117600 56700; -1400 268 ...

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

满秩方程

欠秩方程

超定方程

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

最新文章

热门文章