初始化环境

In [73]:

from IPython.display import display, Math, Latex 
from sympy import *
init_printing()

from helper import comparator_factory, comparator_eval_factory, comparator_method_factory

x,y,z = symbols('x y z')

func_comparator = comparator_factory('使用{}前:','使用后:')
method_comparator = comparator_method_factory('调用{}前:','调用后:')
eval_comparator = comparator_eval_factory('计算前:','计算后:')

矩阵

创建矩阵

从list创建矩阵

如果要在Sympy中创建一个矩阵，可以通过行向量(内嵌的list)构成的list初始化一个Matrix类。

In [74]:

m = Matrix(
    [
        [1, -1], 
        [3, 4], 
        [0, 2]
    ]
)

m

Out[74]:

$$\left[\begin{matrix}1 & -1\\3 & 4\\0 & 2\end{matrix}\right]$$

如果只传入单层的list，那么会创建一个$n \times 1$矩阵。

In [75]:

m = Matrix([1,2,3])

m

Out[75]:

$$\left[\begin{matrix}1\\2\\3\end{matrix}\right]$$

构造常用的矩阵

有一些创建常用矩阵的方法。

单位矩阵

用eye(n)创建单位矩阵。

In [76]:

M = eye(3)

M

Out[76]:

$$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$

零矩阵

通过zero(n)去创建零矩阵。

In [77]:

M = zeros(3,4)

display(M)

$$\left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]$$

Matrix of Ones

通过ones(n)去创建数值全为1的矩阵。

In [78]:

M = ones(3,2)

M

Out[78]:

$$\left[\begin{matrix}1 & 1\\1 & 1\\1 & 1\end{matrix}\right]$$

对角矩阵

传入数个矩阵或者数值去构造对角矩阵。数值会被处理成$1 \times 1$矩阵。

In [79]:

M = diag(-1, ones(2, 2), Matrix([5, 7, 5]))

M

Out[79]:

$$\left[\begin{matrix}-1 & 0 & 0 & 0\\0 & 1 & 1 & 0\\0 & 1 & 1 & 0\\0 & 0 & 0 & 5\\0 & 0 & 0 & 7\\0 & 0 & 0 & 5\end{matrix}\right]$$

基本操作

获得形状

用shape获得矩阵的形状。

In [80]:

M.shape

Out[80]:

$$\left ( 6, \quad 4\right )$$

获取行或者列

使用row()和col()获取单行或者单列。

In [81]:

M.row(0)

Out[81]:

$$\left[\begin{matrix}-1 & 0 & 0 & 0\end{matrix}\right]$$

获取单列。

In [82]:

M.col(-1)

Out[82]:

$$\left[\begin{matrix}0\\0\\0\\5\\7\\5\end{matrix}\right]$$

删除行列

用row_del或col_del删除行或列。这些操作会直接修改矩阵。

例如，我们要删除第一行和最后一列。

In [83]:

M = Matrix([[1, 2, 3], [3, 2, 1]])

display('删除前:',M)

M.row_del(0)
M.col_del(-1)

display('删除后:',M)

'删除前:'

$$\left[\begin{matrix}1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$$

'删除后:'

$$\left[\begin{matrix}3 & 2\end{matrix}\right]$$

矩阵转置

使用T

In [84]:

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

display('原矩阵:',M)
display('转职后矩阵:',M.T)

'原矩阵:'

$$\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\end{matrix}\right]$$

'转职后矩阵:'

$$\left[\begin{matrix}1 & 4\\2 & 5\\3 & 6\end{matrix}\right]$$

插入行列

如果要插入行或者列，用row_insert或col_insert。这些操作并不会直接修改原矩阵而是返回数值被修改的一个新矩阵。

row_insert()接受位置和$1 \times n$矩阵。

In [85]:

M = Matrix([[1, 2, 3], [3, 2, 1]])

method_comparator(M, 'row_insert', 0, Matrix([[6,6,6]]))

调用row_insert()前:

$$\left[\begin{matrix}1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$$

调用后:

$$\left[\begin{matrix}6 & 6 & 6\\1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$$

col_insert()接受位置和$n \times 1$矩阵。

In [86]:

M = Matrix([[1, 2, 3], [3, 2, 1]])

method_comparator(M, 'col_insert', -1,  Matrix([6,6]))

调用col_insert()前:

$$\left[\begin{matrix}1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$$

调用后:

$$\left[\begin{matrix}1 & 2 & 6 & 3\\3 & 2 & 6 & 1\end{matrix}\right]$$

矩阵计算

基本的矩阵计算像加法，减法，乘法，求逆可以通过pythhon的算符+,-,*,**来实现。

加法

In [87]:

M = Matrix([[1, 3], [-2, 3]])
N = Matrix([[0, 3], [0, 7]])

print('M')
display(M)
print('N')
display(N)
print('M+N')
display(M+N)

$$\left[\begin{matrix}1 & 3\\-2 & 3\end{matrix}\right]$$

$$\left[\begin{matrix}0 & 3\\0 & 7\end{matrix}\right]$$

M+N

$$\left[\begin{matrix}1 & 6\\-2 & 10\end{matrix}\right]$$

减法

In [88]:

M = Matrix([[1, 3], [-2, 3]])
N = Matrix([[0, 3], [0, 7]])

print('M')
display(M)
print('N')
display(N)
print('M-N')
display(M-N)

$$\left[\begin{matrix}1 & 3\\-2 & 3\end{matrix}\right]$$

$$\left[\begin{matrix}0 & 3\\0 & 7\end{matrix}\right]$$

M-N

$$\left[\begin{matrix}1 & 0\\-2 & -4\end{matrix}\right]$$

乘法

In [89]:

M = Matrix([[1, 2, 3], [3, 2, 1]])
N = Matrix([0, 1, 1])

print('M')
display(M)
print('N')
display(N)
display('MxN')
display(M*N)

$$\left[\begin{matrix}1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$$

$$\left[\begin{matrix}0\\1\\1\end{matrix}\right]$$

'MxN'

$$\left[\begin{matrix}5\\3\end{matrix}\right]$$

求逆

In [90]:

M = Matrix([[1, 3], [-2, 3]])

print('M')
display(M)
print('inverse of M:')
display(M**(-1))

$$\left[\begin{matrix}1 & 3\\-2 & 3\end{matrix}\right]$$

inverse of M:

$$\left[\begin{matrix}\frac{1}{3} & - \frac{1}{3}\\\frac{2}{9} & \frac{1}{9}\end{matrix}\right]$$

行列式

用det计算矩阵的行列式

In [91]:

M = Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])

method_comparator(M,'det')

调用det()前:

$$\left[\begin{matrix}1 & 0 & 1\\2 & -1 & 3\\4 & 3 & 2\end{matrix}\right]$$

调用后:

$$-1$$

行阶梯型矩阵

用rref()将矩阵处理成行阶梯型矩阵。它返回两个结果，第一个是行阶梯型矩阵，第二个是主元列的index。

In [92]:

M = Matrix([[1, 0, 1, 3], [2, 3, 4, 7], [-1, -3, -3, -4]])

M_rref, pivot_columns = M.rref()

print('原矩阵:')
display(M)

print('阶梯型:')
display(M_rref)

print('主元index:')
display(pivot_columns)

原矩阵:

$$\left[\begin{matrix}1 & 0 & 1 & 3\\2 & 3 & 4 & 7\\-1 & -3 & -3 & -4\end{matrix}\right]$$

阶梯型:

$$\left[\begin{matrix}1 & 0 & 1 & 3\\0 & 1 & \frac{2}{3} & \frac{1}{3}\\0 & 0 & 0 & 0\end{matrix}\right]$$

主元index:

$$\left [ 0, \quad 1\right ]$$

零空间

nullspace()返回一个能够张成矩阵零空间的向量组成的list。

In [93]:

M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])

print('原矩阵M:')
display(M)
print('矩阵M的零空间的张成向量:')
display(M.nullspace())

原矩阵M:

$$\left[\begin{matrix}1 & 2 & 3 & 0 & 0\\4 & 10 & 0 & 0 & 1\end{matrix}\right]$$

矩阵M的零空间的张成向量:

$$\left [ \left[\begin{matrix}-15\\6\\1\\0\\0\end{matrix}\right], \quad \left[\begin{matrix}0\\0\\0\\1\\0\end{matrix}\right], \quad \left[\begin{matrix}1\\- \frac{1}{2}\\0\\0\\1\end{matrix}\right]\right ]$$

列空间

columnspace()返回一个能够张成矩阵列空间的向量组成的list。

In [94]:

M = Matrix([[1, 1, 2], [2 ,1 , 3], [3 , 1, 4]])

print('矩阵M:')
display(M)
print('张成矩阵M的列空间的向量:')
display(M.columnspace())

矩阵M:

$$\left[\begin{matrix}1 & 1 & 2\\2 & 1 & 3\\3 & 1 & 4\end{matrix}\right]$$

张成矩阵M的列空间的向量:

$$\left [ \left[\begin{matrix}1\\2\\3\end{matrix}\right], \quad \left[\begin{matrix}1\\1\\1\end{matrix}\right]\right ]$$

特征值和特征向量

要计算矩阵的特征值，可以使用eigenvals。它返回由特征值:代数重数对组成的dict。

In [95]:

M = Matrix([[3, -2,  4, -2], [5,  3, -3, -2], [5, -2,  2, -2], [5, -2, -3,  3]])

print('矩阵M')
display(M)

print('矩阵M的特征值和它们的代数重数:')
display(M.eigenvals())

矩阵M

$$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$$

矩阵M的特征值和它们的代数重数:

$$\left \{ -2 : 1, \quad 3 : 1, \quad 5 : 2\right \}$$

eigenvectors返回由特征值，几何重数，特征向量列表组成的tuple。

In [96]:

M = Matrix(
    [
        [3, -2,  4, -2], 
        [5,  3, -3, -2], 
        [5, -2,  2, -2], 
        [5, -2, -3,  3]
    ]
)

print('矩阵M')
display(M)

print('特征值，几何重数，特征向量:')
display(M.eigenvects())

矩阵M

$$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$$

特征值，几何重数，特征向量:

$$\left [ \left ( -2, \quad 1, \quad \left [ \left[\begin{matrix}0\\1\\1\\1\end{matrix}\right]\right ]\right ), \quad \left ( 3, \quad 1, \quad \left [ \left[\begin{matrix}1\\1\\1\\1\end{matrix}\right]\right ]\right ), \quad \left ( 5, \quad 2, \quad \left [ \left[\begin{matrix}1\\1\\1\\0\end{matrix}\right], \quad \left[\begin{matrix}0\\-1\\0\\1\end{matrix}\right]\right ]\right )\right ]$$

如果只是想获得特征多项式，用charpoly()和factor()

In [97]:

print('矩阵M')
display(M)

lamda = symbols('lamda')
p = M.charpoly(lamda).factor()

print('矩阵M的特征多项式')
display(p)

矩阵M

$$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$$

矩阵M的特征多项式

$$\left(\lambda - 5\right)^{2} \left(\lambda - 3\right) \left(\lambda + 2\right)$$

对角化

如果要对角化一个矩阵，用diagonalize()。它返回($P$,$D$)组成的tuple，其中$D$是对角矩阵并且满足$M = PDP^{-1}$。

In [98]:

M = Matrix(
    [
        [3, -2,  4, -2], 
        [5,  3, -3, -2], 
        [5, -2,  2, -2], 
        [5, -2, -3,  3]
    ]
)

P, D = M.diagonalize()

print('矩阵M')
display(M)

print('矩阵P')
display(P)

print('矩阵D')
display(D)

display(Latex('$PDP^{-1}$'))
display(P*D*P**-1)

矩阵M

$$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$$

矩阵P

$$\left[\begin{matrix}0 & 1 & 1 & 0\\1 & 1 & 1 & -1\\1 & 1 & 1 & 0\\1 & 1 & 0 & 1\end{matrix}\right]$$

矩阵D

$$\left[\begin{matrix}-2 & 0 & 0 & 0\\0 & 3 & 0 & 0\\0 & 0 & 5 & 0\\0 & 0 & 0 & 5\end{matrix}\right]$$

$PDP^{-1}$

$$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$$

等式求解

基本求解

代数等式求解的主要函数是solveset()。它的语法是solveset(equation, variable=None, domain=S.Complexes)

In [99]:

eq = Eq(x**2-1, 0)

func_comparator(eq, solveset, x)

使用solveset()前:

$$x^{2} - 1 = 0$$

使用后:

$$\left\{-1, 1\right\}$$

如果solveset()中传入一个表达式而非等式，那么它会被假设等于0,。

In [100]:

expr = x**2-1

func_comparator(expr, solveset, x)

使用solveset()前:

$$x^{2} - 1$$

使用后:

$$\left\{-1, 1\right\}$$

求解结果可以是无限集合。

In [101]:

eq = Eq(sin(x) - 1)

func_comparator(eq, solveset, x)

使用solveset()前:

$$\sin{\left (x \right )} - 1 = 0$$

使用后:

$$\left\{2 n \pi + \frac{\pi}{2}\; |\; n \in \mathbb{Z}\right\}$$

求解结果可以是空集。

In [102]:

eq = Eq(exp(x))

func_comparator(eq, solveset, x)

使用solveset()前:

$$e^{x} = 0$$

使用后:

$$\emptyset$$

如果求解失败，会返回一个条件集合。

In [103]:

eq = Eq(cos(x) - x)

func_comparator(eq, solveset, x)

使用solveset()前:

$$- x + \cos{\left (x \right )} = 0$$

使用后:

$$\left\{x\; |\; x \in \mathbb{C} \wedge - x + \cos{\left (x \right )} = 0 \right\}$$

solveset()只会显示每个解一次。

In [104]:

eq = Eq(x**3 - 6*x**2 + 9*x,0)

func_comparator(eq, solveset, x)

使用solveset()前:

$$x^{3} - 6 x^{2} + 9 x = 0$$

使用后:

$$\left\{0, 3\right\}$$

如果要获得多项等式每个解的重数，可以用roots()

In [105]:

eq = Eq(x**3 - 6*x**2 + 9*x,0)

func_comparator(eq, roots, x)

使用roots()前:

$$x^{3} - 6 x^{2} + 9 x = 0$$

使用后:

$$\left \{ 0 : 1, \quad 3 : 2\right \}$$

线性系统

linsolve()用于求解线性系统. 线性系统有不同的定义方法。

等式list

In [106]:

ls = [
    Eq(x + y + z - 1), 
    Eq(x + y + 2*z - 3)
]

func_comparator(ls,linsolve, x, z)

使用linsolve()前:

$$\left [ x + y + z - 1 = 0, \quad x + y + 2 z - 3 = 0\right ]$$

使用后:

$$\left\{\left ( - y - 1, \quad 2\right )\right\}$$

等式list可以简化成表达式list，如果这些表达式都等于0。

In [107]:

ls = [
    x + y + z - 1, 
    x + y + 2*z - 3
]

func_comparator(ls,linsolve, x, z)

使用linsolve()前:

$$\left [ x + y + z - 1, \quad x + y + 2 z - 3\right ]$$

使用后:

$$\left\{\left ( - y - 1, \quad 2\right )\right\}$$

增广矩阵形式

In [108]:

M = Matrix(
    [
        [1, 1, 1, 1], 
        [1, 1, 2, 3]
    ]
)
    
func_comparator(M,linsolve, x, y, z)

使用linsolve()前:

$$\left[\begin{matrix}1 & 1 & 1 & 1\\1 & 1 & 2 & 3\end{matrix}\right]$$

使用后:

$$\left\{\left ( - y - 1, \quad y, \quad 2\right )\right\}$$

$Ax = b$ 形式

In [109]:

M = Matrix(
    [
        [1, 1, 1, 1], 
        [1, 1, 2, 3]
    ]
)
ls = A, b = M[:, :-1], M[:, -1]

func_comparator(ls,linsolve, x, y, z)

使用linsolve()前:

$$\left ( \left[\begin{matrix}1 & 1 & 1\\1 & 1 & 2\end{matrix}\right], \quad \left[\begin{matrix}1\\3\end{matrix}\right]\right )$$

使用后:

$$\left\{\left ( - y - 1, \quad y, \quad 2\right )\right\}$$

非线性系统

solvset()不能解非线性多元系统。遇到这样的需求用solve()代替。

In [110]:

eq_list = [
    Eq(x*y - 1, 0),
    Eq(x-2, 0)
]

func_comparator(eq_list, solve, x, y)

使用solve()前:

$$\left [ x y - 1 = 0, \quad x - 2 = 0\right ]$$

使用后:

$$\left [ \left ( 2, \quad \frac{1}{2}\right )\right ]$$

微分方程

先定义函数符号变量。

In [111]:

f = symbols('f', cls = Function)

然后定义微分方程，例如，$f''(x) - 2f'(x) + f(x) = \sin(x)$

In [112]:

diffeq = Eq(f(x).diff(x, 2) - 2*f(x).diff(x) + f(x), sin(x))

然后用dsolve()求解。

In [113]:

func_comparator(diffeq, dsolve, f(x))

使用dsolve()前:

$$f{\left (x \right )} - 2 \frac{d}{d x} f{\left (x \right )} + \frac{d^{2}}{d x^{2}} f{\left (x \right )} = \sin{\left (x \right )}$$

使用后:

$$f{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{x} + \frac{1}{2} \cos{\left (x \right )}$$

如果无法得到显式的解， dsolve()返回一个隐函数形式。

In [114]:

diffeq = Eq(f(x).diff(x)*(1 - sin(f(x))), 0)

func_comparator(diffeq, dsolve, f(x))

使用dsolve()前:

$$\left(- \sin{\left (f{\left (x \right )} \right )} + 1\right) \frac{d}{d x} f{\left (x \right )} = 0$$

使用后:

$$f{\left (x \right )} + \cos{\left (f{\left (x \right )} \right )} = C_{1}$$

参考资料

Sympy Documentation