sympy-symbolic-math

SymPy — 符号数学

概述

SymPy 是一个用于符号数学的 Python 库，它使用数学符号而非数值近似进行精确计算。它涵盖代数、微积分、方程求解、线性代数、物理和代码生成，全部基于纯 Python 实现，无需外部依赖。

适用场景

• 以符号方式求解方程（代数方程、方程组、微分方程）
• 执行微积分运算（导数、积分、极限、级数展开）
• 化简和操作代数表达式
• 以符号方式处理矩阵（特征值、行列式、分解）
• 将符号表达式转换为快速数值函数（lambdify → NumPy）
• 从数学表达式生成代码（C、Fortran、LaTeX）
• 需要精确结果（例如，sqrt(2) 而非 1.414...）
• 对于数值计算（数组运算、数据的线性代数），请使用 numpy/scipy
• 对于统计建模（回归、假设检验），请使用 statsmodels

前置条件

pip install sympy
# Optional for numerical evaluation:
pip install numpy matplotlib

SymPy 采用纯 Python 实现，没有编译型依赖，可在任何环境中安装。

快速入门

from sympy import symbols, solve, diff, integrate, sqrt, pi

x = symbols('x')

# Solve equation
print(solve(x**2 - 5*x + 6, x))          # [2, 3]

# Derivative
print(diff(x**3 + 2*x, x))                # 3*x**2 + 2

# Integral
print(integrate(x**2, (x, 0, 1)))         # 1/3

# Exact arithmetic
print(sqrt(8))                              # 2*sqrt(2)
print(pi.evalf(30))                        # 3.14159265358979323846264338328

核心 API

1. 符号与表达式

创建符号变量并操作表达式。

from sympy import symbols, Symbol, Rational, S, oo, pi, E, I
from sympy import simplify, expand, factor, collect, cancel, trigsimp

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

# With assumptions (improve simplification)
n = symbols('n', integer=True)
t = symbols('t', positive=True, real=True)
from sympy import sqrt
print(sqrt(t**2))   # t (not Abs(t), because t is positive)

# Exact fractions (avoid floats!)
expr = Rational(1, 3) * x + S(1)/7
print(expr)  # x/3 + 1/7

# Simplification
print(simplify(x**2 + 2*x + 1))          # (x + 1)**2
print(expand((x + 1)**3))                 # x**3 + 3*x**2 + 3*x + 1
print(factor(x**3 - x))                   # x*(x - 1)*(x + 1)
print(collect(x*y + x - 3 + 2*x**2 - z*x**2, x))  # x**2*(2 - z) + x*(y + 1) - 3

2. 微积分

导数、积分、极限和级数。

from sympy import symbols, diff, integrate, limit, series, oo, sin, cos, exp, log

x = symbols('x')

# Derivatives
print(diff(sin(x**2), x))                 # 2*x*cos(x**2)
print(diff(x**4, x, 3))                   # 24*x (third derivative)

# Partial derivatives
x, y = symbols('x y')
f = x**2 * y**3
print(diff(f, x, y))                      # 6*x*y**2

# Integrals
x = symbols('x')
print(integrate(x**2, x))                 # x**3/3 (indefinite)
print(integrate(exp(-x**2), (x, -oo, oo)))  # sqrt(pi) (Gaussian)
print(integrate(x * exp(-x), (x, 0, oo)))  # 1

# Limits
print(limit(sin(x)/x, x, 0))             # 1
print(limit((1 + 1/x)**x, x, oo))        # E

# Taylor series
print(series(exp(x), x, 0, 5))           # 1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)

3. 方程求解

代数方程、超越方程和微分方程。

from sympy import symbols, solve, solveset, Eq, S, linsolve, nonlinsolve, Function, dsolve

x, y = symbols('x y')

# Single equation
print(solve(x**2 - 4, x))                 # [-2, 2]
print(solveset(x**2 - 4, x, S.Reals))     # {-2, 2}

# System of linear equations
print(linsolve([x + y - 5, 2*x - y - 1], x, y))  # {(2, 3)}

# System of nonlinear equations
print(nonlinsolve([x**2 + y - 4, x + y**2 - 4], x, y))

# Differential equation: y'' + y = 0
f = Function('f')
ode = f(x).diff(x, 2) + f(x)
print(dsolve(ode, f(x)))                  # Eq(f(x), C1*sin(x) + C2*cos(x))

# With initial conditions
from sympy import Derivative
ics = {f(0): 1, f(x).diff(x).subs(x, 0): 0}
print(dsolve(ode, f(x), ics=ics))         # Eq(f(x), cos(x))

4. 矩阵与线性代数

符号矩阵运算。

from sympy import Matrix, eye, zeros, ones, diag, symbols

# Create matrices
M = Matrix([[1, 2], [3, 4]])
print(f"Det: {M.det()}")                   # -2
print(f"Inverse:\n{M**-1}")

# Symbolic matrices
a, b = symbols('a b')
M = Matrix([[a, b], [b, a]])
print(f"Eigenvalues: {M.eigenvals()}")     # {a - b: 1, a + b: 1}

# Eigenvectors and diagonalization
eigendata = M.eigenvects()
# [(eigenval, multiplicity, [eigenvectors]), ...]
P, D = M.diagonalize()
print(f"M = P*D*P^-1")

# Solve linear system Ax = b
A = Matrix([[1, 2], [3, 4]])
b = Matrix([5, 6])
x = A.solve(b)
print(f"Solution: {x.T}")

# Matrix calculus
t = symbols('t')
M_t = Matrix([[t, t**2], [1, t]])
print(f"dM/dt:\n{M_t.diff(t)}")

5. 代码生成

将符号表达式转换为快速数值函数或编译型代码。

import numpy as np
from sympy import symbols, lambdify, sin, exp, ccode, fcode, latex

x, y = symbols('x y')
expr = sin(x) * exp(-x**2 / 2)

# lambdify: symbolic → fast NumPy function
f = lambdify(x, expr, 'numpy')
x_vals = np.linspace(-5, 5, 1000)
y_vals = f(x_vals)
print(f"Shape: {y_vals.shape}, Max: {y_vals.max():.4f}")

# Multi-variable lambdify
expr2 = x**2 + y**2
f2 = lambdify((x, y), expr2, 'numpy')
print(f"f(3, 4) = {f2(3, 4)}")            # 25

# C code generation
print(ccode(expr))                          # sin(x)*exp(-1.0/2.0*pow(x, 2))

# Fortran code generation
print(fcode(expr))

# LaTeX output
print(latex(expr))                          # \sin{\left(x \right)} e^{- \frac{x^{2}}{2}}

6. 物理模块

经典力学、向量分析和单位。

from sympy import symbols, cos, sin, Function
from sympy.physics.mechanics import dynamicsymbols, LagrangesMethod, Particle, Point, ReferenceFrame
from sympy.physics.vector import dot, cross

# Vector analysis
N = ReferenceFrame('N')
v1 = 3*N.x + 4*N.y + 0*N.z
v2 = 1*N.x + 0*N.y + 2*N.z
print(f"Dot: {dot(v1, v2)}")               # 3
print(f"Cross: {cross(v1, v2)}")            # 8*N.x - 6*N.y - 4*N.z

# Simple pendulum via Lagrangian mechanics
q = dynamicsymbols('q')       # Generalized coordinate (angle)
m, g, l = symbols('m g l', positive=True)
T = Rational(1, 2) * m * (l * q.diff())**2          # Kinetic energy
V = m * g * l * (1 - cos(q))                         # Potential energy
L = T - V                                            # Lagrangian
print(f"Lagrangian: {L}")

核心概念

精确运算与数值运算

from sympy import Rational, S, sqrt, pi

# WRONG: introduces floating-point error
expr_bad = 0.5 * x           # Float 0.5, loses exactness

# CORRECT: exact symbolic arithmetic
expr_good = Rational(1, 2) * x    # Exact 1/2
expr_good = S(1)/2 * x            # Alternative exact syntax
expr_good = x / 2                  # Also exact

# Numerical evaluation when needed
print(sqrt(2).evalf())         # 1.41421356237310
print(pi.evalf(50))            # 50 digits of precision

求解器选择指南

求解器	适用场景	返回值
solve(eq, x)	通用用途，旧版接口	解的列表
solveset(eq, x, domain)	代数方程（推荐）	集合（可能是无限集）
linsolve(system, vars)	线性方程组	元组的 FiniteSet
nonlinsolve(system, vars)	非线性方程组	元组的 FiniteSet
dsolve(ode, f(x))	常微分方程	等式（Eq）
nsolve(eq, x0)	数值求根	Float 近似值

常见化简函数

函数	作用	示例
simplify()	通用化简（速度较慢，会尝试所有方法）	sin(x)**2 + cos(x)**2 → 1
expand()	展开乘法	(x+1)**2 → x**2+2*x+1
factor()	分解为不可约因式	x**2-1 → (x-1)*(x+1)
collect()	按变量分组	收集 x 的同类项
cancel()	约去分数中的公因式	(x**2-1)/(x-1) → x+1
trigsimp()	化简三角表达式	化简三角表达式时比 simplify 更快
powsimp()	化简幂和指数表达式	合并 x**a * x**b

常见工作流

工作流：符号到数值计算流程

from sympy import symbols, diff, integrate, lambdify, sin, cos
import numpy as np
import matplotlib.pyplot as plt

x = symbols('x')

# 1. Define expression symbolically
f_expr = sin(x) * cos(x)**2

# 2. Symbolic operations
f_prime = diff(f_expr, x)
F_expr = integrate(f_expr, x)
print(f"f(x) = {f_expr}")
print(f"f'(x) = {f_prime}")
print(f"F(x) = {F_expr}")

# 3. Convert to fast numerical functions
f_num = lambdify(x, f_expr, 'numpy')
f_prime_num = lambdify(x, f_prime, 'numpy')
F_num = lambdify(x, F_expr, 'numpy')

# 4. Evaluate and plot
x_vals = np.linspace(0, 2*np.pi, 500)
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].plot(x_vals, f_num(x_vals)); axes[0].set_title('f(x)')
axes[1].plot(x_vals, f_prime_num(x_vals)); axes[1].set_title("f'(x)")
axes[2].plot(x_vals, F_num(x_vals)); axes[2].set_title('F(x)')
plt.tight_layout()
plt.savefig('symbolic_pipeline.png', dpi=150)
print("Saved symbolic_pipeline.png")

工作流：求解并验证

from sympy import symbols, solve, simplify, Eq, sqrt

x = symbols('x')

# 1. Define equation
equation = x**3 - 6*x**2 + 11*x - 6

# 2. Solve symbolically
solutions = solve(equation, x)
print(f"Solutions: {solutions}")            # [1, 2, 3]

# 3. Verify each solution
for sol in solutions:
    result = simplify(equation.subs(x, sol))
    assert result == 0, f"Solution {sol} failed!"
    print(f"  x={sol}: f(x) = {result} ✓")

# 4. Factor the polynomial
from sympy import factor
print(f"Factored: {factor(equation)}")      # (x - 1)*(x - 2)*(x - 3)

工作流：ODE 系统分析

1. 使用 Function 和 dsolve() 定义 ODE
2. 以符号方式求解；使用 ics={} 参数应用初始条件
3. 使用 lambdify() 将解转换为数值函数
4. 使用 matplotlib 绘制解的轨迹

关键参数

参数	函数	默认值	选项	作用
domain	solveset()	S.Complexes	S.Reals、S.Integers	限定解的定义域
force	simplify()	False	True/False	强力化简
n	diff(expr, x, n)	1	1–∞	导数阶数
精度	evalf(n)	15	1–1000+	数值精度位数
后端	lambdify()	"math"	"numpy"、"scipy"、"mpmath"	数值计算后端
rational	nsimplify()	True	True/False	求精确有理数近似值

最佳实践

1. 始终使用 Rational() 或 S() 表示分数：0.5 * x 会引入浮点数，破坏精确计算。请使用 Rational(1, 2) * x 或 S(1)/2 * x。

2. 为符号添加假设：symbols('x', positive=True) 可启用诸如 sqrt(x**2) → x 的化简。没有假设时，SymPy 必须处理一般复数情况。

3. 使用 lambdify 进行数值求值，而不是 subs().evalf()：在循环中使用 subs/evalf，速度比单次 lambdify 调用慢 100–1000 倍。

   # Slow: [expr.subs(x, v).evalf() for v in values]
   # Fast: f = lambdify(x, expr, 'numpy'); f(np.array(values))
   

4. 反模式——默认使用 simplify()：simplify() 会尝试许多策略，因此速度较慢。如果知道目标形式，请使用特定函数（factor、expand、trigsimp）。

5. 代数方程优先使用 solveset，而不是 solve：solveset 会返回合适的数学集合，并能更好地处理边界情况。solve 属于旧版接口，但在一般情况下仍然有用。

6. 反模式——数值求解已足够时仍进行符号求解：对于没有闭式解的方程，请使用 nsolve(eq, x0) 进行数值求根，而不是等待 solve 失败。

7. 在 Jupyter 中使用 init_printing()，以获得易读的输出：from sympy import init_printing; init_printing() 可在笔记本中启用 LaTeX 渲染。

常用方案

方案：生成 LaTeX 文档

from sympy import symbols, Integral, Eq, latex, sqrt, pi

x = symbols('x')
integral = Integral(x**2 * sqrt(1 - x**2), (x, 0, 1))
result = integral.doit()

print(f"$$ {latex(integral)} = {latex(result)} $$")
# $$ \int\limits_{0}^{1} x^{2} \sqrt{1 - x^{2}}\, dx = \frac{\pi}{16} $$

方案：参数化 ODE 解

from sympy import symbols, Function, dsolve, Eq, exp, lambdify
import numpy as np

x = symbols('x')
k, A = symbols('k A', positive=True)
f = Function('f')

# Solve y' = -ky with y(0) = A
ode = Eq(f(x).diff(x), -k * f(x))
solution = dsolve(ode, f(x), ics={f(0): A})
print(f"Solution: {solution}")              # f(x) = A*exp(-k*x)

# Evaluate for specific parameters
f_num = lambdify((x, k, A), solution.rhs, 'numpy')
x_vals = np.linspace(0, 5, 100)
y_vals = f_num(x_vals, k=0.5, A=10)
print(f"y(5) = {y_vals[-1]:.4f}")

方案：符号矩阵分解

from sympy import Matrix, symbols, pprint

a, b, c, d = symbols('a b c d')
M = Matrix([[a, b], [c, d]])

# Characteristic polynomial
lam = symbols('lambda')
char_poly = M.charpoly(lam)
print(f"Characteristic polynomial: {char_poly.as_expr()}")

# Eigenvalues (symbolic)
eigenvals = M.eigenvals()
print(f"Eigenvalues: {eigenvals}")

# Determinant and trace
print(f"det(M) = {M.det()}")               # a*d - b*c
print(f"tr(M) = {M.trace()}")              # a + d

故障排除

问题	原因	解决方案
NameError: name 'x' is not defined	未创建符号	使用前通过 x = symbols('x') 定义
结果中出现意外的浮点数	使用了 0.5，而不是 Rational(1,2)	使用 Rational() 或 S() 表示精确分数
simplify() 速度很慢	尝试对复杂表达式应用所有策略	使用特定函数：factor()、expand()、trigsimp()
solve() 返回空列表	不存在闭式解	使用 nsolve(eq, x0) 进行数值近似
sqrt(x**2) 返回 sqrt(x**2)，而不是 x	未对 x 添加假设	使用 x = symbols('x', positive=True) 定义
lambdify 结果错误	表达式包含 SymPy 专用函数	指定后端：lambdify(x, expr, 'numpy') 或 'scipy'
dsolve 中出现 NotImplementedError	不支持此类 ODE	改用数值 ODE 求解器（scipy odeint）

相关技能

• matplotlib-scientific-plotting — 将符号结果转换为 lambdify 形式后进行绘图
• statsmodels-statistical-modeling — 统计推断；需要 p 值而非精确代数运算时使用
• matlab-scientific-computing — 用于数值（而非符号）计算的 MATLAB 替代方案

参考资料

• SymPy 文档 — 官方 API 参考和教程
• SymPy 教程 — 面向新用户的入门介绍
• SymPy GitHub — 源代码、示例和问题
• Meurer 等人（2017）《SymPy：Python 中的符号计算》— PeerJ Computer Science