6.1 Matrix 与 ImmutableMatrix
SymPy 的矩阵体系有两个核心类:Matrix 和 ImmutableMatrix。它们共享同一套运算接口,区别只有一个——可变性。Matrix 是可变对象,你可以原地修改它的元素;ImmutableMatrix 是不可变对象,任何“修改”都返回新矩阵,原对象不变。这个区别在单次脚本里看起来无关紧要,但在符号计算中影响深远:可变对象不可哈希。
先看创建和可变性:
PYTHON
import sympy as sp
A = sp.Matrix([[1, 2], [3, 4]])
print(A)
print("形状:", A.shape)
# 可变:原地修改元素
A[0, 0] = 99
print("修改后:", A)
# 不可变
B = sp.ImmutableMatrix([[1, 2], [3, 4]])
try:
B[0, 0] = 99
except TypeError as e:
print("ImmutableMatrix 修改失败:", e)
输出
Matrix([[1, 2], [3, 4]])
形状: (2, 2)
修改后: Matrix([[99, 2], [3, 4]])
ImmutableMatrix 修改失败: 'Cannot set values of 'sympy.matrices.immutable.ImmutableDenseMatrix'>
Matrix 接受嵌套列表、扁平列表加维度、行列数加生成函数等多种输入。对数学工作者来说最自然的是嵌套列表:外层是行,内层是列。如果传入扁平列表,必须显式指定 (rows, cols),否则 SymPy 会按方阵或列向量猜测——猜错时不报错,只给你形状不对的矩阵。
可变性的实际后果是哈希:
PYTHON
A = sp.Matrix([[1, 2], [3, 4]])
B = sp.ImmutableMatrix([[1, 2], [3, 4]])
try:
d = {A: "mutable"}
except TypeError as e:
print("Matrix 作字典键失败:", e)
d = {B: "immutable"}
print("ImmutableMatrix 作字典键成功:", d[B])
输出
Matrix 作字典键失败: unhashable type: 'MutableDenseMatrix'
ImmutableMatrix 作字典键成功: immutable
这在符号计算中很重要:假设系统、表达式缓存、subs 的替换字典,底层都依赖哈希。如果你要把矩阵嵌入更复杂的符号结构(矩阵作为表达式系数、矩阵作为函数参数缓存),ImmutableMatrix 是唯一安全的选择。本书大部分示例用 Matrix,因为交互式探索中可变性更方便;工程代码中 ImmutableMatrix 更常见。
转换是显式的,没有隐式转换,两套类型不能直接混用运算符:
PYTHON
A = sp.Matrix([[1, 2], [3, 4]])
B = A.as_immutable()
C = B.as_mutable()
print(type(A).__name__, type(B).__name__, type(C).__name__)
输出
MutableDenseMatrix ImmutableDenseMatrix MutableDenseMatrix
这不是 SymPy 的缺陷,是可变/不可变二分法的固有代价:编译器无法在运行时替你决定“这个操作应该原地改还是返回新对象”。
还有一个容易忽略的点:Matrix 的元素可以是符号表达式,但矩阵本身不是表达式。Matrix.subs() 逐元素替换,返回新矩阵;Eq 对矩阵的比较是结构比较,不是数学相等。这些细节在 6.6 节会反复出现。
6.2 行列式、逆、秩与零空间
这一节的核心函数是 inv、rank、null、trace、transpose、ctranspose。它们覆盖矩阵最基本的“体检指标”。SymPy 对每个函数都提供符号和数值两条路径,行为差异值得逐一说明。
求逆 inv
Matrix.inv() 默认走符号精确路径:先算行列式,非零则用伴随矩阵或高斯消元给出精确逆。对数值矩阵,精确路径返回有理数矩阵,不是浮点数矩阵:
PYTHON
A = sp.Matrix([[1, 2], [3, 4]])
A_inv = A.inv()
print("精确逆:", A_inv)
print("元素类型:", [type(x).__name__ for x in A_inv])
print("数值逆:", A_inv.evalf())
print("验证 A * A_inv:", A * A_inv)
输出
精确逆: Matrix([[-2, 1], [3/2, -1/2]])
元素类型: ['Integer', 'One', 'Rational', 'Rational']
数值逆: Matrix([[-2.00000000000000, 1.00000000000000], [1.50000000000000, -0.500000000000000]])
验证 A * A_inv: Matrix([[1, 0], [0, 1]])
evalf() 把精确逆转成浮点逆,改变了返回类型。对教学演示没问题,但对“判断矩阵是否奇异”这类问题,evalf 之后的逆矩阵不再携带“精确奇异”的信息。本站的 qikjik_inv 在返回前调了 .evalf(),读者需要知道这个封装在精度语义上做了简化。
inv 的 method 参数可选 'GE'、'LU'、'ADJ'、'CH'、'LDL'。不同方法可能给出形式上不同但数学等价的逆矩阵;'CH' 要求对称正定,'LDL' 要求厄米特,不满足时抛错。这些是数学约束,不是实现约束,封装没有替你检查。
秩 rank
Matrix.rank() 对符号和数值矩阵都走精确路径,返回整数:
PYTHON
A = sp.Matrix([[1, 2], [2, 4]]) # 第二行是第一行的 2 倍
print("rank:", A.rank())
# 符号矩阵
a = sp.symbols('a')
M = sp.Matrix([[1, 2], [2, a]])
print("符号 rank:", M.rank()) # 对一般 a 为 2
print("a=4 时 det:", M.subs(a, 4).det()) # 0,秩降为 1
输出
rank: 1
符号 rank: 2
a=4 时 det: 0
本站的 qikjik_rank 做了一个重要分支:数值矩阵走 SVD + 容差,符号矩阵走 Matrix.rank()。动机是对标 MATLAB 的 rank:浮点矩阵的“秩”在数值意义上依赖容差,精确行化简会把舍入误差放大成“满秩”。SVD 路径下小于容差的奇异值被视为零。读者需要知道:同一个浮点矩阵,qikjik_rank 和 Matrix.rank() 可能给出不同结果;符号矩阵没有这个问题。
零空间 null
Matrix.nullspace() 返回列向量列表,基是有理基(精确):
PYTHON
A = sp.Matrix([[1, 2], [3, 6]])
ns = A.nullspace()
print("零空间基:", ns)
print("验证 A * v = 0:", A * ns[0])
输出
零空间基: [Matrix([-2, 1])]
验证 A * v = 0: Matrix([0, 0])
本站的 qikjik_null 提供两条路径:tol='rational' 时调 nullspace() 返回有理基;tol 为浮点数时走 SVD 返回标准正交基。有理基和正交基张成同一个零空间,但基向量不同:有理基元素是有理数,适合符号推导;正交基向量模长为 1、相互正交,适合数值投影和最小二乘。两者不能混用。
迹 trace
Matrix.trace() 返回对角线元素之和,对非方阵取 min(rows, cols) 个对角元素:
PYTHON
A = sp.Matrix([[1, 2], [3, 4]])
print("trace:", A.trace())
print("非方阵 trace:", sp.Matrix([[1, 2, 3], [4, 5, 6]]).trace())
输出
trace: 5
sympy.matrices.common.NonSquareMatrixError
本站的 qikjik_trace 额外处理标量输入:对标量返回标量本身。SymPy 原生 trace 不接受标量。封装价值在于兼容 MATLAB 习惯。
转置与共轭转置
.T 是普通转置,.H 是共轭转置:
PYTHON
A = sp.Matrix([[1 + 2*sp.I, 3], [4, 5 - sp.I]])
print("A.T:", A.T)
print("A.H:", A.H)
输出
A.T: Matrix([[1 + 2*I, 4], [3, 5 - I]])
A.H: Matrix([[1 - 2*I, 4], [3, 5 + I]])
对实数矩阵两者相同,对复数矩阵 .H 先转置再取共轭。这里要提醒一个常见混淆:MATLAB 的 A' 是共轭转置,A.' 是普通转置;SymPy 的 .H 对应 A',.T 对应 A.'。本站函数名 ctranspose 指向共轭转置,transpose 指向普通转置,命名清晰,但从 MATLAB 迁移时容易把 ' 的语义带错。
6.3 特征值与特征向量:符号计算的独特优势
特征值问题是符号计算真正发光的地方。数值库给你浮点特征值,符号库给你精确代数表达式,包括带参数的特征值。这个差异在工程中经常被低估:当你需要知道“特征值随某个参数如何变化”时,数值特征值只能给你离散采样点,符号特征值给你一个函数。
SymPy 的接口是 eigenvals() 和 eigenvects()。前者返回字典,键是特征值,值是代数重数:
PYTHON
A = sp.Matrix([[1, 2], [3, 4]])
eigs = A.eigenvals()
print("特征值字典:", eigs)
# 重根示例
B = sp.Matrix([[2, 0], [0, 2]])
print("重根:", B.eigenvals()) # {2: 2}
输出
特征值字典: {5/2 - sqrt(33)/2: 1, 5/2 + sqrt(33)/2: 1}
重根: {2: 2}
eigenvects() 返回三元组列表 (特征值, 重数, 特征向量列表):
PYTHON
A = sp.Matrix([[1, 2], [3, 4]])
for val, mult, vecs in A.eigenvects():
print(f"λ={val}, 代数重数={mult}, 特征向量={vecs}")
输出
λ=5/2 - sqrt(33)/2, 代数重数=1, 特征向量=[Matrix([...])]
λ=5/2 + sqrt(33)/2, 代数重数=1, 特征向量=[Matrix([...])]
注意 eigenvects 返回的是几何重数对应的特征向量,不是代数重数对应的广义特征向量。如果矩阵不可对角化,几何重数小于代数重数,eigenvects 给出的向量数少于矩阵阶数。这是 6.4 节 Jordan 标准型的伏笔。
本站的 qikjik_eig 做了三层封装:n=1 只返回特征值(展开重数,返回列表);n=2 返回右特征向量矩阵和对角阵;n=3 额外返回左特征向量矩阵。n=2 和 n=3 走 Matrix.diagonalize(),要求矩阵可对角化,不可对角化时抛 ValueError。这意味着 qikjik_eig(A, 2) 对亏损矩阵直接失败,而不是像 MATLAB 的 eig 那样返回广义特征向量或 Jordan 形式。
符号特征值的独特优势在于参数化:
PYTHON
a = sp.symbols('a', real=True)
M = sp.Matrix([[a, 1], [1, a]])
print("含参特征值:", M.eigenvals())
输出
含参特征值: {a - 1: 1, a + 1: 1}
你可以直接分析“当 a 取何值时特征值变为重根”“特征值的符号如何随 a 变化”。这是数值特征值做不到的:数值库对每个 a 给一组数,你要自己拟合趋势。第 15 章符号-数值混合优化会回到这个主题:Hessian 的符号特征值告诉你临界点的稳定性,数值特征值只能告诉你“在这个点上稳定或不稳定”,不能告诉你“在参数空间的哪个区域稳定”。
另一个细节:n=3 的左特征向量实现是 matrix.T.diagonalize()[0].T。左特征向量满足 wTA = λwT,等价于 ATw = λw,所以左特征向量矩阵是 AT 的右特征向量矩阵的转置。这个推导是对的,实现依赖 AT 也可对角化——如果 A 可对角化,AT 也可对角化,依赖是安全的。
6.4 对角化与 Jordan 标准型
对角化的数学条件是:矩阵有 n 个线性无关的特征向量,等价于几何重数等于代数重数对每个特征值成立。Matrix.diagonalize() 返回 (P, D),满足 A = P D P-1:
PYTHON
A = sp.Matrix([[1, 2], [3, 4]])
P, D = A.diagonalize()
print("P:", P)
print("D:", D)
print("验证 P*D*P^-1:", sp.simplify(P * D * P.inv()))
输出
P: Matrix([[-sqrt(33)/6 - 1/2, -1/2 + sqrt(33)/6], [1, 1]])
D: Matrix([[5/2 - sqrt(33)/2, 0], [0, 5/2 + sqrt(33)/2]])
验证 P*D*P^-1: Matrix([[1, 2], [3, 4]])
不可对角化时抛 ValueError:
PYTHON
J = sp.Matrix([[1, 1], [0, 1]]) # Jordan 块
try:
J.diagonalize()
except ValueError as e:
print("不可对角化:", e)
输出
不可对角化: Matrix is not diagonalizable
rref(简化行阶梯形)是对角化的底层工具之一,也是判断秩、解线性方程组、求零空间的通用工具。Matrix.rref() 返回 (R, pivots):
PYTHON
A = sp.Matrix([[1, 2, 1, 1], [2, 4, 2, 2], [1, 1, 1, 3]])
R, pivots = A.rref()
print("R:", R)
print("pivots:", pivots)
输出
R: Matrix([[1, 0, 1, 5], [0, 1, 0, -2], [0, 0, 0, 0]])
pivots: (0, 1)
本站的 qikjik_rref 只返回 R,丢弃了 pivots。这个取舍对“看结果”够用,对“分析结构”不够:主元列索引告诉你哪些变量是自由变量,这是解线性方程组的关键信息。读者如果需要完整信息,直接用原生 rref()。
compan(伴随矩阵)是一个容易被低估的工具。给定首一多项式 p(x) = xn + a1xn-1 + ... + an,伴随矩阵的特征值恰好是 p 的根:
PYTHON
import sympy as sp
import numpy as np
def qikjik_compan(u):
# 处理多项式表达式输入
if isinstance(u, sp.Expr) and u.free_symbols:
x = tuple(u.free_symbols)[0]
poly = sp.Poly(u, x)
raw_coeffs = poly.all_coeffs()
# 归一化处理:使首项系数为1
leading_coeff = raw_coeffs[0]
coeffs = [c / leading_coeff for c in raw_coeffs[1:]]
# 处理列表或矩阵输入
else:
if isinstance(u, sp.Matrix):
coeffs = list(u.flat())
else:
coeffs = list(u)
n = len(coeffs)
matrix = np.zeros((n, n), dtype=object)
# 填充次对角线
for i in range(n - 1):
matrix[i + 1, i] = 1
# 填充第一行负系数
matrix[0, :] = [-c for c in coeffs]
return sp.Matrix(matrix)
x = sp.symbols('x')
p = sp.Poly(x**2 - 3*x + 2)
C = sp.Matrix.companion(p) if hasattr(sp.Matrix, 'companion') else None
# 或用本站 qikjik_compan
C = qikjik_compan([1, -3, 2])
print("伴随矩阵:", C)
print("特征值:", C.eigenvals()) # {1: 1, 2: 1},即 p 的根
输出
伴随矩阵: Matrix([[-1, 3, -2], [1, 0, 0], [0, 1, 0]])
特征值: {-(3*sqrt(321)/2 + 83/2)**(1/3)/3 - 10/(3*(3*sqrt(321)/2 + 83/2)**(1/3)) - 1/3: ...
这把“求多项式根”转化为“求矩阵特征值”,在数值线性代数中是标准技巧(companion matrix + QR 算法)。本站 qikjik_compan 对表达式输入先用 Poly.all_coeffs() 提取系数再除以首项系数归一化;对列表输入不做归一化,直接取列表作为系数。这个行为对标 MATLAB 的 compan(u):要求 u 是首一多项式的系数,不自动归一化。读者需要自己保证输入是首一的。
Jordan 标准型在 SymPy 中没有稳定的顶层函数。实际做法是:先用 eigenvects() 判断是否可对角化;不可对角化时,用 rref 求广义特征向量链。本书不展开完整算法,因为它的符号计算复杂度高,且数值上病态。第 11 章会讨论这个性能问题。
6.5 符号矩阵的局限与性能考量
符号矩阵的第一个局限是表达式膨胀。行列式的精确展开是指数级的:
PYTHON
n = 5
M = sp.Matrix(n, n, lambda i, j: sp.Symbol(f'a{i}{j}'))
det = M.det()
print("5x5 符号行列式项数:", sp.count_ops(det))
输出
5x5 符号行列式项数: 599
3×3 行列式展开有 6 项,4×4 有 24 项,5×5 有 120 项。inv 的伴随矩阵法需要算 n² 个代数余子式,每个余子式是一个 (n-1)×(n-1) 行列式。对 5×5 矩阵,符号求逆的中间表达式规模可能达到数万项。这不是 SymPy 的缺陷,是符号代数的固有复杂度。
第二个局限是化简的不可判定性。符号矩阵的运算结果往往需要 simplify 才能变成可读形式,但 simplify 不保证成功,也不保证最快。第 3 章已经讨论过这个问题,矩阵场景只是把它放大了:矩阵元素的化简需要逐个进行,而元素之间的代数关系(比如行列式为零)可能跨越多个元素,单元素化简看不到全局结构。
cond(条件数)和 norm(范数)是数值工具,不是符号工具:
PYTHON
import sympy as sp
import numpy as np
def qikjik_cond(A, p=2):
if isinstance(A, sp.Matrix):
np_matrix = np.array(A, dtype=float)
else:
np_matrix = np.array(A, dtype=float)
# 计算条件数
if p == 'inf' or p == np.inf:
p = np.inf
elif p == 'fro':
p = 'fro'
return np.linalg.cond(np_matrix, p=p)
A = sp.Matrix([[4.1, 2.8], [9.7, 6.6]])
print("2-范数条件数:", sp.N(qikjik_cond(A))) # 走浮点路径
print("Frobenius 范数:", A.norm('fro'))
输出
2-范数条件数: 1622.99938385651
Frobenius 范数: 12.7397017233529
本站的 qikjik_cond 和 qikjik_norm 都走 NumPy/SciPy 路径,返回浮点数。这是合理的:条件数的定义依赖奇异值,奇异值的符号表达式对一般矩阵不可行。读者需要知道:符号矩阵的 cond 和 norm 是数值近似,不是精确值;对符号矩阵调用这两个函数,内部会先转成浮点数组,精度损失不可避免。如果你需要“符号条件数”,只能对特殊结构的矩阵(对角、三角)手工推导。
norm 的 p 参数在 Matrix.norm() 中支持 1、2、oo、'fro'。本站封装透传这些选项,但对向量和矩阵的语义没有区分:norm([1, 2, 3]) 在 SymPy 中会把列表转成 3×1 矩阵,然后算矩阵范数,结果和向量 2-范数一致,但这不是显式的向量范数接口。
6.6 实战:符号化线性系统分析
一个完整的符号化线性系统分析案例:给定含参数 a 的线性方程组,分析解的存在性、唯一性和解的结构。
方程组:
x + 2y + z = 1
2x + 4y + 2z = 2
x + y + a*z = 3
第一步,构造增广矩阵并做 rref:
PYTHON
import sympy as sp
a, x, y, z = sp.symbols('a x y z')
A = sp.Matrix([[1, 2, 1, 1], [2, 4, 2, 2], [1, 1, a, 3]])
R, pivots = A.rref()
print("R:", R)
print("pivots:", pivots)
输出
R: Matrix([[1, 0, 2*a - 1, 5], [0, 1, 1 - a, -2], [0, 0, 0, 0]])
pivots: (0, 1)
第二行被消成零行(因为第二行是第一行的 2 倍),第一、三行给出主元。pivots 告诉你主元列,自由变量是剩余列对应的变量。这个分析是精确的,a 保留为符号。
第二步,用 linsolve 求解:
PYTHON
M = sp.Matrix([[1, 2, 1], [2, 4, 2], [1, 1, a]])
b = sp.Matrix([1, 2, 3])
sol = sp.linsolve((M, b), x, y, z)
print(sol)
输出
{(z*(1 - 2*a) + 5, z*(a - 1) - 2, z)}
返回的 FiniteSet 会包含 a 的表达式,以及在 a 取特定值(使矩阵奇异)时的特殊情形。linsolve 对参数化线性系统的处理是符号计算的优势:它给出“对所有 a 成立”的解结构,而不是对每个 a 数值求解再拼凑。
第三步,分析系数矩阵的秩和行列式:
PYTHON
print("det:", M.det()) # 0(第二行是第一行的2倍)
print("rank:", M.rank()) # 对大多数 a 为 2
输出
det: 0
rank: 2
行列式为零说明矩阵奇异,解不唯一或不存在。秩为 2 说明零空间维数为 1,解空间是一条直线(如果相容)或空集(如果不相容)。这些结论是符号推导的直接结果,不需要数值采样。
第四步,特征值分析(如果矩阵是方阵且你关心稳定性):
PYTHON
sp.pprint(M.eigenvals())
输出
⎧ ________________ ________________ ⎫
⎪ ╱ 2 ╱ 2 ⎪
⎨ a ╲╱ a - 10⋅a + 37 5 a ╲╱ a - 10⋅a + 37 5 ⎬
⎪0: 1, ─ - ─────────────────── + ─: 1, ─ + ─────────────────── + ─: 1⎪
⎩ 2 2 2 2 2 2 ⎭
对含参数的矩阵,特征值表达式告诉你“参数取何值时出现重根”“特征值的符号如何变化”。这是第 15 章优化项目的基础:Hessian 的特征值符号决定临界点的稳定性。
完整链条:rref 看结构 → linsolve 求解 → det / rank 判断存在性 → eigenvals 分析稳定性。每一步都是符号的,参数 a 贯穿始终。数值方法做不到这一点:你必须对每个 a 采样,然后拟合趋势,而且无法保证没有漏掉临界值。
6.7 本站对照
以下函数是 MATLAB 风格接口的 SymPy 参考实现。多数矩阵分解走 SciPy/NumPy 数值路径,不是符号分解。
矩阵分解
-
chol(A):Cholesky 分解,要求对称正定。
PYTHON
# 原生:SymPy 有 Matrix.cholesky(),限制多
# 本站:走 scipy.linalg.cholesky
qikjik_chol([[4, 2], [2, 3]]) # 返回 L
见 https://www.qikjik.com/math_func_examples/?p=chol(A)
-
lu(A):LU 分解,返回 L、U、P 满足 A = P'LU。
PYTHON
# 原生:Matrix.LUdecomposition() 不返回 MATLAB 形式的 P
# 本站:走 scipy.linalg.lu
qikjik_lu([[1, 2], [3, 4]])
见 https://www.qikjik.com/math_func_examples/?p=lu(A)
-
qr(A):QR 分解,返回上三角 R。
PYTHON
# 原生:Matrix.QRdecomposition() 存在,符号路径慢
# 本站:走 scipy.linalg.qr
qikjik_qr([[1, 2], [3, 4]])
见 https://www.qikjik.com/math_func_examples/?p=qr(A)
-
svd(A):奇异值分解,返回 U、S、V。
PYTHON
# 原生:SymPy 无符号 SVD
# 本站:走 numpy.linalg.svd
qikjik_svd([[1, 2], [3, 4]])
见 https://www.qikjik.com/math_func_examples/?p=svds(A)
-
pinv(A):Moore-Penrose 伪逆。
PYTHON
# 原生:Matrix.pinv() 存在(符号精确,慢)
# 本站:走 numpy.linalg.pinv(数值)
qikjik_pinv([[1, 2], [3, 4]])
见 https://www.qikjik.com/math_func_examples/?p=pinv(A)
这些分解在 SymPy 中没有完整的符号实现(Matrix.cholesky() 存在但限制多,Matrix.LUdecomposition() 存在但不返回 MATLAB 形式的置换矩阵)。本站实现走数值路径,返回浮点矩阵。读者需要知道:符号分解和数值分解是不同的问题;符号分解给出精确表达式但可能膨胀,数值分解给出浮点结果但稳定。
矩阵函数
-
expm(X):矩阵指数。
PYTHON
# 原生:SymPy 无 Matrix.exp() 符号实现(需先对角化)
# 本站:走 scipy.linalg.expm
qikjik_expm([[0, 1], [-1, 0]])
见 https://www.qikjik.com/math_func_examples/?p=expm(X)
-
logm(A):主矩阵对数。
PYTHON
# 原生:无
# 本站:走 scipy.linalg.logm
qikjik_logm([[1, 0], [0, 2]])
见 https://www.qikjik.com/math_func_examples/?p=logm(A)
-
sqrtm(A):主矩阵平方根。
PYTHON
# 原生:无
# 本站:走 scipy.linalg.sqrtm
qikjik_sqrtm([[4, 0], [0, 9]])
见 https://www.qikjik.com/math_func_examples/?p=sqrtm(A)
-
mpower(A,B):矩阵幂 A^B。
PYTHON
# 原生:Matrix**n 对整数幂有符号实现
# 本站:mpower 支持整数幂,非整数幂走数值
qikjik_mpower([[1, 1], [0, 1]], 3)
见 https://www.qikjik.com/math_func_examples/?p=mpower(A,B)
矩阵函数在 SymPy 中的符号支持有限:exp 对矩阵没有原生符号实现,需要先对角化再逐元素取指数。本站实现走数值路径(SciPy 的 expm、logm、sqrtm),返回浮点矩阵。对符号矩阵,这些函数会先转浮点,精度损失不可避免。
矩阵运算
-
kron(A,B):Kronecker 张量积。
PYTHON
# 原生
sp.Matrix([[1, 2]]).kronecker_product(sp.Matrix([[3], [4]]))
# 本站
qikjik_kron([[1, 2]], [[3], [4]])
见 https://www.qikjik.com/math_func_examples/?p=kron(A,B)
-
dot(A,B):点积。
PYTHON
# 原生
sp.Matrix([1, 2, 3]).dot(sp.Matrix([4, 5, 6]))
# 本站
qikjik_dot([1, 2, 3], [4, 5, 6])
见 https://www.qikjik.com/math_func_examples/?p=dot(A,B)
-
cross(A,B):叉积。
PYTHON
# 原生
sp.Matrix([1, 0, 0]).cross(sp.Matrix([0, 1, 0]))
# 本站
qikjik_cross([1, 0, 0], [0, 1, 0])
见 https://www.qikjik.com/math_func_examples/?p=cross(A,B)
kron 在 SymPy 中有原生 Matrix.kronecker_product(),符号矩阵也支持。dot 和 cross 对符号向量有原生实现(Matrix.dot()、Matrix.cross()),本站封装主要是输入格式兼容。
矩阵结构判断
bandwidth(A, type):上下带宽。
isbanded / isdiag / ishermitian / issymmetric / istril / istriu:结构判断。
PYTHON
A = sp.Matrix([[1, 2], [0, 3]])
print("上三角:", A.is_upper) # True
print("对角:", sp.Matrix([[1, 0], [0, 2]]).is_diagonal) # True
print("对称:", sp.Matrix([[1, 2], [2, 3]]).is_symmetric) # True
这些判断对符号矩阵有效,因为只看元素位置是否为零。ishermitian 对符号矩阵需要比较 A 和 A.H 的元素是否恒等,SymPy 的 is_hermitian 属性会做这个比较,但对含参数的矩阵可能返回 None(无法判定)。读者需要知道:结构判断的返回值可能是 True、False 或 None,None 表示“无法判定”,不是“否”。
向量微积分
hessian(f, v):Hessian 矩阵。
jacobian(f, v):Jacobian 矩阵。
laplacian(f, v):Laplacian 算子。
gradient(F):梯度。
divergence(X,Y,Z,U,V,W):散度。
curl(V):旋度。
potential(V,X):向量场势函数。
surfnorm(function,varlist,point):曲面法向量。
PYTHON
x, y = sp.symbols('x y')
f = x**2 + x*y + y**2
# 原生
print("hessian:", sp.hessian(f, (x, y)))
print("jacobian:", sp.Matrix([f]).jacobian([x, y]))
# 本站
qikjik_hessian(f, [x, y])
qikjik_jacobian([f], [x, y])
这些函数是第 14 章单摆项目和后续物理推导的基础。hessian 和 jacobian 在 SymPy 中有原生实现(Matrix.jacobian(),以及 sympy.vector 模块的 hessian),但本站封装提供了 MATLAB 风格的调用接口。gradient、divergence、curl 在 sympy.vector 模块中有完整实现,但需要先定义坐标系和向量场对象,比 MATLAB 的扁平调用更繁琐。本站封装的价值在于降低入门门槛,代价是隐藏了 sympy.vector 的坐标系抽象。