SymPy中的向量和矩阵

SymPy 中的向量和矩阵

Vectors and Matrice in SymPy

准备工作：

from sympy import *
init_printing()
a,b,c,d,e,f,g,h,i = symbols('a:i') # neat shorthand for multiple symbols!

带符号的矩阵

1、一般用矩阵类 Matrix([[],[]]) 生成可以改变的矩阵(mutable matrix)

A = Matrix([[a,b,c],[d,e,f]])
A
A[i,j]  # a(i,j), i=0,1,..., j=0,1,...
A[:,0], A[1,:] # 第一列向量j=0，第二行向量i=1
B = A[:,0:2]
B[0,0] = i

$B=\left[\begin{array}{cc}i& b\\ d& e\end{array}\right]$

2、还可以简便操作，提供行列数和数列

B = Matrix(2,3,[1,2,3,4,5,6])

$B=\left[\begin{array}{ccc}1& 2& 3\\ 4& 5& 6\end{array}\right]$

3、通过二元函数来定义矩阵

>>> B=Matrix(3,4,lambda i,j: 1-(i+j)%2)

$B=\left[\begin{array}{cccc}1& 0& 1& 0\\ 0& 1& 0& 1\\ 1& 0& 1& 0\end{array}\right]$

4、用命令 ImmutableMatrix() 生成不可以更改的矩阵(immutable matrix)

C = ImmutableMatrix([[a, b,],[c, d]])
C[0, 0] = i # 报错，不可更改矩阵元素的值

5、如何得到矩阵的行列数

获取矩阵属性值：行数、列数 M.rows, M.cols, M.shape

>>> A.rows, A.cols, A.shape #(行数=2, 列数=3, 行列数=(2, 3))

6、转置矩阵(Transpose of a matrix)

A.T

$A=\left[\begin{array}{cc}a& d\\ b& e\\ c& f\end{array}\right]$

7、行向量 row vector 和列向量 column vector

单独一行或一列就是行向量或列向量。
没有特别定义行向量和列向量，同样用Matrix命令来定义。

Matrix([[a, b]]) #定义行向量 $[ab]$
Matrix([[c], [d]]) 或者 Matrix([c,d]) #定义列向量 $\left[\begin{array}{l}c\\ d\end{array}\right]$

RVector = Matrix([[a,b,c]) # 行向量
CVector1 = Matrix([[c],[d],[e]]) # 列向量1
CVector2 = Matrix([c,d,e]) # 列向量2

$RVector=[a,b,c]$,
$CVector=\left[\begin{array}{c}a\\ b\\ c\end{array}\right]$

通过指定行列数来获取行、列向量，更加方便快捷

>>> V=Matrix(1,3,[x,y,z])  #行向量，一行，V=[x  y  z]
>>> V=Matrix(3,1,[x,y,z])  #列向量，一列，V=[x  y  z].T

>>> M
Matrix([
[a, b, c],
[d, e, f]])
>>> M[0:M.rows*M.cols]  #组合成一行
[a, b, c, d, e, f]
>>> M[M.cols:M.cols*M.rows]  #取第二行
[d, e, f]

矩阵运算(Matrix Operations)

1、线代 (linear algebra) 核心是向量的线性运算

包括加+, 减-, 乘*(内积或点积 inner or dot product, v_1.T*v_2, v_1.dot(v_2)), 外积(outer product, v_1*v_2.T), 除/。

v_1, v_2 = Matrix([c,d]), Matrix([e,f])
M_1= a*v_1+b*v_2 # 线性加法
M_2 = v_1.T*v_2  # 点积, 结果是仅仅含有一个数的矩阵
v = v_1.dot(v_2) # 点积得到的结果是一个数值
u = v_1*v_2.T    # 外积得一个矩阵

$M_1=\left[\begin{array}{c}ac+be\\ ad+bf\end{array}\right]$ $M_2=[ce+df]v=ce+dfu=\left[\begin{array}{cc}ce& cf\\ de& df\end{array}\right]$

2、求偏导 M.diff(x)

V = Matrix([x,y,1]) #列向量，齐次坐标
M = A*V        #矩阵乘积
D = M.diff(x)  #对x求偏导

$M=\left[\begin{array}{c}ax+by+c\\ dx+ey+f\end{array}\right]$ $D=\left[\begin{array}{c}a\\ d\end{array}\right]$

3、逆矩阵 M**-1

>>> C=Matrix([[1,2],[3,4]])
>>> C**-1  #逆矩阵inverse
Matrix([
[-2,     1],
[3/2, -1/2]])

4、方阵求行列式的值 M.det()

>>> C.det()  #方阵的行列式求值
-2

5、添删矩阵的某行或某列

在某些情况下，需要对矩阵进行加上一行或者加上一列的操作，在这里就可以使用 row_insert 或者
col_insert 函数：第一个参数表示插入的位置，第二个参数就是相应的行向量或者列向量。

而在删除上就很简单了，直接使用 row_del 或者 col_del 即可.

>>> A.row_insert(1,Matrix([[10,10]]))  #第2行上插入
Matrix([
[ 1,  0],
[10, 10],
[ 0,  1]])
>>> A.col_insert(0,Matrix([3,5]))  #第1列上插入
Matrix([
[3, 1, 0],
[5, 0, 1]])
>>> A.row_del(0)  #删除第0+1行
>>> A
Matrix([[0, 1]])
>>> A.col_del(1)  #删除第1+1=2列
>>> A
Matrix([[0]])

6、叉积cross product

叉积采用 M.cross(A) 的方式实现

v_1, v_2 = Matrix([1,2,3]), Matrix([4,5,6])
v_3 = v_1.cross(v_2)
v_1.dot(v_3)  #0 因为 v_1和v_2都 垂直 v_3
v_2.dot(v_3)  #0

7、函数对矩阵元素的操作 M.applyfunc(myfunc)

f = lambda x: 2*x  #定义单变量函数
eye(3).applyfunc(f) #作用到矩阵的每个元素上

8、替换矩阵中的变量 M.subs(x,4)

M = eys(3)*x
M.subs(x,5) #将x用5代替
M.subs(x,y) #将x替换成y

矩阵属性(Matrix Properties)

SymPy提供了一个判断矩阵是否是带符号的矩阵的方法 .is_symbolic()

A.is_symbolic()  # 返回 True
C = Matrix[[1,2],[3,4]]
C.is_symbolic()  # 返回 False

创建特殊矩阵

• 单位矩阵 eye(n,m)
  在矩阵论中，最常见的就是单位矩阵了，而单位矩阵只与一个参数有关，那就是矩阵的大小, Sympy中用函数 eye(n,m) 创建 $n\times m$的单位矩阵。

>>> eye(4)
Matrix([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])
>>> eye(2,3)
Matrix([
[1, 0, 0],
[0, 1, 0]])

• 零矩阵和全一矩阵 zeros(n,m), ones(n,m)
  另外还有全零 zeros(n,m) 和 全一 ones(n,m) 矩阵，意思就是矩阵中的所有值全部是 零和一。

>>> ones(3,4)
Matrix([
[1, 1, 1, 1],
[1, 1, 1, 1],
[1, 1, 1, 1]])
>>> zeros(3,2)
Matrix([
[0, 0],
[0, 0],
[0, 0]])

• 对角矩阵 diag(n,m,k,...)
  而对角矩阵也可以使用 diag(n,m,k,...) 轻松获得：

>>> diag(1,2,3,4)
Matrix([
[1, 0, 0, 0],
[0, 2, 0, 0],
[0, 0, 3, 0],
[0, 0, 0, 4]])

• 矩阵特征值、特征向量与对角化
  在对角化方面，同样可以使用 特征值eigenvals()，特征向量eigenvects(), 对角化diagonalize() 函数

>>> A=sympy.Matrix([[1,1],[0,2]])
>>> A
Matrix([
[1, 1],
[0, 2]])
>>> A.eigenvals()   # 特征值
{1: 1, 2: 1}
>>> A.eigenvects()  # 特征向量
[(1, 1, [Matrix([
[1],
[0]])]), (2, 1, [Matrix([
[1],
[1]])])]
>>> A.diagonalize()  # 对角化
(Matrix([
[1, 1],
[0, 1]]), Matrix([
[1, 0],
[0, 2]]))

>>> A=Matrix([[1,2],[2,2]])
>>> A.eigenvals()
{3/2 - sqrt(17)/2: 1, 3/2 + sqrt(17)/2: 1}
>>> A.eigenvects()
[(3/2 - sqrt(17)/2, 1, [Matrix([
[-2/(-1/2 + sqrt(17)/2)],
[                     1]])]), (3/2 + sqrt(17)/2, 1, [Matrix([
[-2/(-sqrt(17)/2 - 1/2)],
[                     1]])])]
>>> A.diagonalize()
(Matrix([
[-sqrt(17)/4 - 1/4, -1/4 + sqrt(17)/4],
[                1,                 1]]), Matrix([
[3/2 - sqrt(17)/2,                0],
[               0, 3/2 + sqrt(17)/2]]))

在 eigenvals() 返回的结果中，第一个表示特征值，第二个表示该特征值的重数。在特征向量 eigenvects()
中，第一个表示特征值，第二个表示特征值的重数，第三个表示特征向量。
在对角化 diagonalize() 中，第一个矩阵表示 $P$，第二个矩阵表示 $D,A=P\times D\times P^{-1}$。

• 矩阵方程求解
  在矩阵中，最常见的还是多元一次方程的解。如果要求 $Ax=b$ 的解，可以有以下方案:

>>> A
Matrix([
[1, 1],
[0, 2]])
>>> b = sympy.Matrix([3,5])
>>> b
Matrix([
[3],
[5]])
>>> A**-1*b
Matrix([
[1/2],
[5/2]])
>>> sympy.linsolve((A,b))
FiniteSet((1/2, 5/2))
>>> sympy.linsolve((A,b),[x,y])
FiniteSet((1/2, 5/2))

创建符号矩阵的过程

def make_Aij(m, n, a='a') :
    from sympy import Symbol, Matrix # just in case they aren't already loaded
    A = zeros(m, n)
    for i in range(0, m) :
        for j in range(0, n) :
            s = a+'_'+str(i)+str(j)
            exec s + "= Symbol('" + s + "')" # go look up what "exec" does!
            exec "A[i, j] = " + s
    return A
C = make_Aij(2, 3, 'c')

创建矩阵 $C=\left[\begin{array}{ccc}{c}_{00}& c_{01}& c_{02}\\ c_{10}& c_{11}& c_{12}\end{array}\right]$

这里的 exec 函数用来立刻生成符号（字符串）。