sympy

SymPy - 符号数学

用于符号数学的 Python 库，完全以 Python 提供计算机代数系统（CAS）功能。

适用场景

• 符号表达式和代数运算
• 微积分（导数、积分、极限、级数展开）
• 求解方程（代数方程、超越方程、微分方程）
• 表达式的化简和变换
• 矩阵运算和线性代数（符号计算）
• 特殊数学函数
• 数学证明和验证
• 代码生成（C、Fortran、LaTeX）
• 物理计算（力学、量子力学）
• 数论和离散数学
• 逻辑和布尔代数
• 几何代数

参考文档

官方文档：https://docs.sympy.org/
搜索模式：sympy.symbols、sympy.diff、sympy.integrate、sympy.solve、sympy.simplify

核心原则

使用 SymPy 处理

任务	模块	示例
符号	symbols	x, y = symbols('x y')
导数	diff	diff(x**2, x)
积分	integrate	integrate(x**2, x)
求解方程	solve	solve(x**2 - 4, x)
化简	simplify	simplify(expr)
极限	limit	limit(sin(x)/x, x, 0)
级数展开	series	series(exp(x), x, 0, 5)
矩阵	Matrix	Matrix([[1, 2], [3, 4]])

不要用于

• 数值计算（使用 NumPy、SciPy）
• 快速数值计算（符号计算速度较慢）
• 机器学习（使用 PyTorch、TensorFlow）
• 统计分析（使用 SciPy、statsmodels）
• 大规模数值模拟（使用 NumPy、Numba）

快速参考

安装

# pip
pip install sympy

# conda
conda install sympy

# With all optional dependencies
pip install sympy[all]

标准导入

# Core imports
import sympy as sp
from sympy import symbols, Symbol, Function
from sympy import diff, integrate, limit, series

# Common operations
from sympy import simplify, expand, factor, collect
from sympy import solve, solveset, dsolve

# Special functions
from sympy import sin, cos, exp, log, sqrt, Abs
from sympy import pi, E, I, oo  # Constants

# Matrices
from sympy import Matrix, eye, zeros, ones

# Printing
from sympy import init_printing, pprint, latex
init_printing()  # Pretty printing in notebook

基本模式 - 符号表达式

from sympy import symbols, simplify, expand

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

# Create expression
expr = (x + y)**2

# Expand
expanded = expand(expr)
print(f"Expanded: {expanded}")  # x**2 + 2*x*y + y**2

# Factor back
from sympy import factor
factored = factor(expanded)
print(f"Factored: {factored}")  # (x + y)**2

# Substitute values
result = expr.subs([(x, 2), (y, 3)])
print(f"Result: {result}")  # 25

基本模式 - 微积分

from sympy import symbols, diff, integrate, limit
from sympy import sin, cos, exp

x = symbols('x')

# Derivative
f = x**3 + 2*x**2 - x + 1
df = diff(f, x)
print(f"f'(x) = {df}")  # 3*x**2 + 4*x - 1

# Integral
integral = integrate(sin(x), x)
print(f"∫sin(x)dx = {integral}")  # -cos(x)

# Definite integral
definite = integrate(x**2, (x, 0, 1))
print(f"∫₀¹ x²dx = {definite}")  # 1/3

# Limit
lim = limit(sin(x)/x, x, 0)
print(f"lim(sin(x)/x) as x→0 = {lim}")  # 1

关键规则

✅ 应该做

• 使用 symbols() 定义变量 - 正确地定义符号变量
• 化简表达式 - 使用 simplify() 清理结果
• 使用有理数 - 使用 Rational(1, 2) 代替 0.5
• 检查假设条件 - 必要时为符号设置假设条件
• 使用 subs() 进行代换 - 用值替换符号
• 美化打印结果 - 使用 pprint() 或 init_printing()
• 通过数值计算验证结果 - 转换为浮点数进行检查
• 使用合适的函数 - 选择正确的求解函数
• 求解前先因式分解 - 先化简方程
• 使用 lambdify 提高速度 - 转换为 NumPy 函数

❌ 不要做

• 不谨慎地混用符号和数值 - 明确指定类型
• 在符号计算中使用 Python 浮点数 - 使用 Rational 或 Integer
• 忘记定义符号 - 使用前必须先声明
• 忽略符号与数值的区别 - 知道何时使用各自的计算方式
• 使用 == 求解方程 - 使用 solve() 或 Eq()
• 盲目执行开销大的运算 - 某些积分很难求解
• 假设会自动化简 - 通常需要显式调用 simplify()
• 用符号计算处理大型数值任务 - 速度太慢
• 忘记假设条件 - 这可能影响结果（正数、实数等）
• 过度依赖 solve() - 使用 solveset() 以获得更好的处理效果

反模式（绝对不要）

from sympy import symbols, solve, simplify, integrate, Rational
import sympy as sp

# ❌ BAD: Using float instead of Rational
x = symbols('x')
expr = x + 0.5  # Float!
# Result may not be exact

# ✅ GOOD: Use Rational for exact arithmetic
expr = x + Rational(1, 2)
# Exact representation

# ❌ BAD: Not defining symbols
result = y**2 + 2*y + 1  # NameError: y not defined

# ✅ GOOD: Define symbols first
y = symbols('y')
result = y**2 + 2*y + 1

# ❌ BAD: Using Python == for equations
solve(x**2 == 4)  # Wrong! Returns boolean

# ✅ GOOD: Use solve() properly
solve(x**2 - 4, x)  # Returns [-2, 2]

# Or use Eq()
from sympy import Eq
solve(Eq(x**2, 4), x)

# ❌ BAD: Not simplifying
expr = (x + 1)**2 - (x**2 + 2*x + 1)
print(expr)  # Messy, not simplified to 0

# ✅ GOOD: Simplify
result = simplify(expr)
print(result)  # 0

# ❌ BAD: Using symbolic for numerical loops
for i in range(1000000):
    result = sp.sin(sp.pi * i / 180)  # Very slow!

# ✅ GOOD: Lambdify for numerical work
import numpy as np
x = symbols('x')
f_sym = sp.sin(x)
f_num = sp.lambdify(x, f_sym, 'numpy')
results = f_num(np.linspace(0, np.pi, 1000000))  # Fast!

# ❌ BAD: Ignoring assumptions
x = symbols('x')
sqrt(x**2)  # Returns sqrt(x**2), not x (could be negative)

# ✅ GOOD: Set assumptions
x = symbols('x', positive=True)
simplify(sqrt(x**2))  # Returns x

符号和表达式

创建符号

from sympy import symbols, Symbol

# Single symbol
x = Symbol('x')

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

# With assumptions
x = symbols('x', real=True)
y = symbols('y', positive=True)
n = symbols('n', integer=True)
theta = symbols('theta', real=True)

# Complex symbols
z = symbols('z', complex=True)

# Functions
from sympy import Function
f = Function('f')
g = Function('g')

# Indexed symbols
from sympy import IndexedBase, Idx
A = IndexedBase('A')
i, j = symbols('i j', integer=True)
element = A[i, j]

print(f"Assumptions for x: {x.assumptions0}")

构建表达式

from sympy import symbols, sin, cos, exp, log, sqrt
from sympy import pi, E, I, oo

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

# Arithmetic
expr1 = x + 2*y - 3*z
expr2 = x**2 + y**2
expr3 = x*y / z

# Functions
expr4 = sin(x) + cos(y)
expr5 = exp(x**2)
expr6 = log(x + 1)

# Constants
expr7 = pi * x
expr8 = E**x
expr9 = I * x  # Imaginary unit

# Complex expressions
expr10 = (x + y)**3 / (sqrt(x**2 + y**2))
expr11 = sin(x)**2 + cos(x)**2

# Accessing parts of expressions
print(f"Numerator: {expr3.as_numer_denom()[0]}")
print(f"Denominator: {expr3.as_numer_denom()[1]}")

# Expression info
print(f"Free symbols: {expr10.free_symbols}")
print(f"Is polynomial: {expr1.is_polynomial()}")

代换

from sympy import symbols, sin, cos, pi

x, y = symbols('x y')

# Basic substitution
expr = x**2 + 2*x + 1
result = expr.subs(x, 3)
print(f"Result: {result}")  # 16

# Multiple substitutions
expr = x**2 + y**2
result = expr.subs([(x, 1), (y, 2)])
print(f"Result: {result}")  # 5

# Substitute expression
expr = sin(x) + cos(x)
result = expr.subs(x, pi/4)
print(f"Result: {result}")  # sqrt(2)

# Sequential substitution
expr = x + y
temp = expr.subs(x, y)
result = temp.subs(y, 1)
print(f"Result: {result}")  # 2

# Substitute and simplify
from sympy import simplify
expr = sin(x)**2 + cos(x)**2
result = simplify(expr.subs(x, pi/3))
print(f"Result: {result}")  # 1

微积分

导数

from sympy import symbols, diff, sin, cos, exp, log
from sympy import pprint

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

# First derivative
f = x**3 + 2*x**2 - x + 1
df = diff(f, x)
print(f"f'(x) = {df}")  # 3*x**2 + 4*x - 1

# Higher order derivatives
d2f = diff(f, x, 2)  # Second derivative
d3f = diff(f, x, 3)  # Third derivative
print(f"f''(x) = {d2f}")
print(f"f'''(x) = {d3f}")

# Partial derivatives
f = x**2*y + y**3*z
df_dx = diff(f, x)
df_dy = diff(f, y)
df_dz = diff(f, z)
print(f"∂f/∂x = {df_dx}")  # 2*x*y
print(f"∂f/∂y = {df_dy}")  # x**2 + 3*y**2*z

# Mixed partial derivatives
d2f_dxdy = diff(f, x, y)
d2f_dydz = diff(f, y, z)
print(f"∂²f/∂x∂y = {d2f_dxdy}")  # 2*x

# Derivative of special functions
g = sin(x**2) * exp(-x)
dg = diff(g, x)
print(f"g'(x) = {dg}")

# Chain rule automatically applied
h = sin(cos(x))
dh = diff(h, x)
print(f"h'(x) = {dh}")  # -sin(x)*cos(cos(x))

积分

from sympy import symbols, integrate, sin, cos, exp, log, sqrt
from sympy import pi, oo

x, y = symbols('x y')

# Indefinite integrals
f = x**2
F = integrate(f, x)
print(f"∫x²dx = {F}")  # x**3/3

# Multiple integrals
f = x*y
F = integrate(f, x)
print(f"∫xy dx = {F}")  # x**2*y/2

# Definite integrals
result = integrate(x**2, (x, 0, 1))
print(f"∫₀¹ x²dx = {result}")  # 1/3

# Improper integrals
result = integrate(exp(-x), (x, 0, oo))
print(f"∫₀^∞ e^(-x)dx = {result}")  # 1

# Trigonometric integrals
result = integrate(sin(x)**2, x)
print(f"∫sin²(x)dx = {result}")  # x/2 - sin(2*x)/4

# Integration by parts (automatic)
result = integrate(x*exp(x), x)
print(f"∫x·e^x dx = {result}")

# Multiple integration
f = x*y
result = integrate(f, (x, 0, 1), (y, 0, 2))
print(f"∫∫xy dx dy = {result}")  # 1

# Difficult integrals
from sympy import erf
result = integrate(exp(-x**2), (x, -oo, oo))
print(f"∫₋∞^∞ e^(-x²)dx = {result}")  # sqrt(pi)

极限

from sympy import symbols, limit, sin, cos, exp, log
from sympy import oo, pi

x = symbols('x')

# Basic limits
lim = limit(sin(x)/x, x, 0)
print(f"lim(sin(x)/x) as x→0 = {lim}")  # 1

# Limits at infinity
lim = limit((1 + 1/x)**x, x, oo)
print(f"lim((1+1/x)^x) as x→∞ = {lim}")  # E

# One-sided limits
lim_left = limit(1/x, x, 0, '-')
lim_right = limit(1/x, x, 0, '+')
print(f"Left limit: {lim_left}")   # -oo
print(f"Right limit: {lim_right}")  # oo

# Limits with L'Hôpital's rule (automatic)
lim = limit((exp(x) - 1)/x, x, 0)
print(f"lim((e^x - 1)/x) as x→0 = {lim}")  # 1

# Multivariable limits
y = symbols('y')
lim = limit(limit(x*y/(x**2 + y**2), x, 0), y, 0)
print(f"Double limit: {lim}")

# Limit of sequences
n = symbols('n', integer=True, positive=True)
lim = limit((1 + 1/n)**n, n, oo)
print(f"lim((1+1/n)^n) as n→∞ = {lim}")  # E

级数展开

from sympy import symbols, series, sin, cos, exp, log
from sympy import O  # Big-O notation

x = symbols('x')

# Taylor series
s = series(exp(x), x, 0, 5)
print(f"exp(x) ≈ {s}")
# 1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)

# Remove O term
s_no_o = s.removeO()
print(f"Without O term: {s_no_o}")

# Series around different points
s = series(log(x), x, 1, 5)
print(f"log(x) around x=1: {s}")

# Trigonometric series
s = series(sin(x), x, 0, 7)
print(f"sin(x) ≈ {s}")
# x - x**3/6 + x**5/120 + O(x**7)

# Laurent series (negative powers)
s = series(1/(x**2 + x), x, 0, 3)
print(f"1/(x²+x) ≈ {s}")

# Series substitution
s1 = series(exp(x), x, 0, 5)
s2 = series(sin(x), x, 0, 5)
s_composed = series(s1.removeO().subs(x, s2.removeO()), x, 0, 5)
print(f"exp(sin(x)) ≈ {s_composed}")

# Multivariable series
y = symbols('y')
s = series(exp(x*y), x, 0, 3)
print(f"exp(xy) in x: {s}")

求解方程

代数方程

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

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

# Simple equation
solutions = solve(x**2 - 4, x)
print(f"x² = 4: {solutions}")  # [-2, 2]

# Using Eq() for clarity
eq = Eq(x**2, 4)
solutions = solve(eq, x)
print(f"Solutions: {solutions}")

# Multiple solutions
eq = x**3 - 6*x**2 + 11*x - 6
solutions = solve(eq, x)
print(f"x³ - 6x² + 11x - 6 = 0: {solutions}")  # [1, 2, 3]

# Quadratic formula (symbolic)
a, b, c = symbols('a b c')
eq = a*x**2 + b*x + c
solutions = solve(eq, x)
print(f"ax² + bx + c = 0:")
for sol in solutions:
    print(f"  x = {sol}")

# Rational equations
eq = 1/x + 1/(x+1) - 1/2
solutions = solve(eq, x)
print(f"Solutions: {solutions}")

# Equations with radicals
eq = sqrt(x) + sqrt(x - 1) - 2
solutions = solve(eq, x)
print(f"Solutions: {solutions}")

方程组

from sympy import symbols, solve, Eq

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

# Linear system
eq1 = Eq(x + y, 5)
eq2 = Eq(x - y, 1)
solution = solve([eq1, eq2], [x, y])
print(f"Linear system: {solution}")  # {x: 3, y: 2}

# Nonlinear system
eq1 = x**2 + y**2 - 4
eq2 = x - y - 1
solutions = solve([eq1, eq2], [x, y])
print(f"Nonlinear system: {solutions}")

# Three equations, three unknowns
eq1 = x + y + z - 6
eq2 = 2*x - y + z - 2
eq3 = x + 2*y - z - 2
solution = solve([eq1, eq2, eq3], [x, y, z])
print(f"3x3 system: {solution}")

# Parametric solutions
eq1 = x + 2*y + z - 1
eq2 = 2*x + y + 2*z - 2
# Underdetermined system
solution = solve([eq1, eq2], [x, y, z])
print(f"Parametric solution: {solution}")

微分方程

from sympy import symbols, Function, dsolve, Eq
from sympy import diff, sin, cos, exp

x = symbols('x')
f = Function('f')

# First-order ODE
# f'(x) = f(x)
eq = Eq(diff(f(x), x), f(x))
solution = dsolve(eq, f(x))
print(f"f'(x) = f(x): {solution}")  # f(x) = C1*exp(x)

# Second-order ODE
# f''(x) + f(x) = 0
eq = Eq(diff(f(x), x, 2) + f(x), 0)
solution = dsolve(eq, f(x))
print(f"f''(x) + f(x) = 0: {solution}")
# f(x) = C1*sin(x) + C2*cos(x)

# ODE with initial conditions
eq = Eq(diff(f(x), x), f(x))
ics = {f(0): 1}  # Initial condition: f(0) = 1
solution = dsolve(eq, f(x), ics=ics)
print(f"Solution with IC: {solution}")  # f(x) = exp(x)

# Non-homogeneous ODE
# f''(x) + 4*f(x) = sin(x)
eq = Eq(diff(f(x), x, 2) + 4*f(x), sin(x))
solution = dsolve(eq, f(x))
print(f"Non-homogeneous: {solution}")

# System of ODEs
g = Function('g')
eq1 = Eq(diff(f(x), x), g(x))
eq2 = Eq(diff(g(x), x), -f(x))
solution = dsolve([eq1, eq2], [f(x), g(x)])
print(f"System of ODEs: {solution}")

超越方程和数值求解

from sympy import symbols, solve, sin, cos, exp, log
from sympy import nsolve, lambdify
import numpy as np

x = symbols('x')

# Transcendental equation
eq = sin(x) - x/2
solutions = solve(eq, x)
print(f"sin(x) = x/2 solutions: {solutions}")

# Some equations can't be solved symbolically
# Use nsolve for numerical solution
eq = exp(x) + x
# nsolve requires initial guess
try:
    solution = nsolve(eq, -1)  # Start near -1
    print(f"e^x + x = 0: x ≈ {solution}")
except:
    print("Could not find solution")

# Trigonometric equation
eq = sin(x) + cos(x) - 1
solutions = solve(eq, x)
print(f"sin(x) + cos(x) = 1: {solutions[:3]}")  # First 3

# Logarithmic equation
eq = log(x) - 1/x
solutions = solve(eq, x)
print(f"log(x) = 1/x: {solutions}")

化简和运算

化简

from sympy import symbols, simplify, sin, cos, exp, log, sqrt
from sympy import pi

x, y = symbols('x y')

# Basic simplification
expr = (x**2 + 2*x + 1) / (x + 1)
simplified = simplify(expr)
print(f"Simplified: {simplified}")  # x + 1

# Trigonometric simplification
expr = sin(x)**2 + cos(x)**2
simplified = simplify(expr)
print(f"sin²x + cos²x = {simplified}")  # 1

# Exponential simplification
expr = exp(x) * exp(y)
simplified = simplify(expr)
print(f"e^x · e^y = {simplified}")  # exp(x + y)

# Radical simplification
expr = sqrt(x**2)
x_pos = symbols('x', positive=True)
simplified = simplify(expr.subs(x, x_pos))
print(f"√(x²) = {simplified}")  # x (if x positive)

# Complex simplification
expr = (x + y)**3 - (x**3 + 3*x**2*y + 3*x*y**2 + y**3)
simplified = simplify(expr)
print(f"Result: {simplified}")  # 0

# Logarithmic simplification
expr = log(x) + log(y)
simplified = simplify(expr)
print(f"log(x) + log(y) = {simplified}")  # log(x*y)

展开和因式分解

from sympy import symbols, expand, factor, collect
from sympy import sin, cos

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

# Expand polynomial
expr = (x + y)**3
expanded = expand(expr)
print(f"(x+y)³ = {expanded}")
# x**3 + 3*x**2*y + 3*x*y**2 + y**3

# Expand trigonometric
expr = sin(x + y)
expanded = expand(expr, trig=True)
print(f"sin(x+y) = {expanded}")
# sin(x)*cos(y) + sin(y)*cos(x)

# Factor polynomial
expr = x**2 + 2*x + 1
factored = factor(expr)
print(f"x² + 2x + 1 = {factored}")  # (x + 1)**2

# Factor complex expressions
expr = x**4 - 1
factored = factor(expr)
print(f"x⁴ - 1 = {factored}")
# (x - 1)*(x + 1)*(x**2 + 1)

# Collect terms
expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
collected = collect(expr, x)
print(f"Collected: {collected}")
# x**3 + x**2*(2 - z) + x*(y + 1) - 3

# Expand and factor
expr = (x + 1)*(x + 2)*(x + 3)
expanded = expand(expr)
refactored = factor(expanded)
print(f"Expanded: {expanded}")
print(f"Factored: {refactored}")

重写表达式

from sympy import symbols, sin, cos, exp, log, tan
from sympy import sinh, cosh

x = symbols('x')

# Rewrite in terms of different functions
expr = tan(x)
rewritten = expr.rewrite(sin)
print(f"tan(x) in terms of sin: {rewritten}")
# sin(x)/cos(x)

# Exponential form
expr = sin(x)
rewritten = expr.rewrite(exp)
print(f"sin(x) in exp form: {rewritten}")
# -I*(exp(I*x) - exp(-I*x))/2

# Hyperbolic to exponential
expr = sinh(x)
rewritten = expr.rewrite(exp)
print(f"sinh(x) = {rewritten}")
# (exp(x) - exp(-x))/2

# Logarithm base change
expr = log(x, 10)  # log base 10
rewritten = expr.rewrite(log)  # Natural log
print(f"log₁₀(x) = {rewritten}")

# Trigonometric identities
expr = sin(2*x)
rewritten = expr.rewrite(sin, cos)
print(f"sin(2x) = {rewritten}")  # 2*sin(x)*cos(x)

部分分式

from sympy import symbols, apart, together
from sympy import cancel, factor

x = symbols('x')

# Partial fraction decomposition
expr = (x**2 + 2*x + 1) / (x**3 + x**2)
partial = apart(expr, x)
print(f"Partial fractions: {partial}")
# 1/x + 2/x**2 - 1/(x + 1)

# Recombine fractions
expr = 1/x + 2/x**2 - 1/(x + 1)
combined = together(expr)
print(f"Combined: {combined}")

# Cancel common factors
expr = (x**2 - 1) / (x - 1)
cancelled = cancel(expr)
print(f"Cancelled: {cancelled}")  # x + 1

# More complex example
expr = (4*x**3 + 21*x**2 + 10*x + 12) / (x**4 + 5*x**3 + 5*x**2 + 4*x)
partial = apart(expr, x)
print(f"Complex partial fractions: {partial}")

线性代数

矩阵

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

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

# Create matrices
A = Matrix([[1, 2], [3, 4]])
B = Matrix([[x, y], [z, x+y]])

print(f"Matrix A:\n{A}")

# Special matrices
I = eye(3)  # Identity
Z = zeros(2, 3)  # Zeros
O = ones(3, 2)  # Ones
D = diag(1, 2, 3)  # Diagonal

print(f"Identity:\n{I}")
print(f"Diagonal:\n{D}")

# Matrix operations
C = A + A
D = A * A  # Matrix multiplication
E = 2 * A

print(f"A + A:\n{C}")
print(f"A * A:\n{D}")

# Transpose
At = A.T
print(f"A transpose:\n{At}")

# Symbolic matrices
M = Matrix([[x, y], [y, z]])
print(f"Symbolic matrix:\n{M}")

矩阵性质

from sympy import Matrix, symbols
from sympy import sqrt

# Create matrix
A = Matrix([[1, 2, 3],
            [4, 5, 6],
            [7, 8, 9]])

# Determinant
det_A = A.det()
print(f"det(A) = {det_A}")  # 0 (singular)

# Rank
rank_A = A.rank()
print(f"rank(A) = {rank_A}")  # 2

# Trace
trace_A = A.trace()
print(f"trace(A) = {trace_A}")  # 15

# Inverse (if exists)
B = Matrix([[1, 2], [3, 4]])
B_inv = B.inv()
print(f"B inverse:\n{B_inv}")

# Verify: B * B^(-1) = I
product = B * B_inv
print(f"B * B⁻¹:\n{product}")

# Eigenvalues and eigenvectors
eigenvals = B.eigenvals()
print(f"Eigenvalues: {eigenvals}")

eigenvects = B.eigenvects()
print(f"Eigenvectors:")
for val, mult, vects in eigenvects:
    print(f"  λ = {val}, multiplicity = {mult}")
    for v in vects:
        print(f"    v = {v}")

矩阵分解

from sympy import Matrix, eye
from sympy import QR, LU

# QR decomposition
A = Matrix([[1, 2], [3, 4], [5, 6]])
Q, R = A.QRdecomposition()
print(f"Q:\n{Q}")
print(f"R:\n{R}")

# Verify: A = Q*R
print(f"Q*R:\n{Q*R}")

# LU decomposition
B = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
L, U, perm = B.LUdecomposition()
print(f"L:\n{L}")
print(f"U:\n{U}")

# Jordan form
from sympy import Matrix
M = Matrix([[2, 1], [0, 2]])
P, J = M.jordan_form()
print(f"Jordan form J:\n{J}")
print(f"Transform P:\n{P}")

# Diagonalization (if possible)
A = Matrix([[1, 2], [2, 1]])
try:
    P, D = A.diagonalize()
    print(f"Diagonal form D:\n{D}")
    print(f"Transform P:\n{P}")
    # Verify: A = P*D*P^(-1)
    print(f"P*D*P⁻¹:\n{P*D*P.inv()}")
except:
    print("Matrix not diagonalizable")

求解线性方程组

from sympy import Matrix, symbols

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

# System: Ax = b
A = Matrix([[3, 2, -1],
            [2, -2, 4],
            [-1, 0.5, -1]])

b = Matrix([1, -2, 0])

# Solve using matrix inverse
x_sol = A.inv() * b
print(f"Solution using inverse:\n{x_sol}")

# Solve using solve_linear_system
from sympy import solve_linear_system
# Augmented matrix
M = Matrix([[3, 2, -1, 1],
            [2, -2, 4, -2],
            [-1, 0.5, -1, 0]])

solution = solve_linear_system(M, x, y, z)
print(f"Solution: {solution}")

# Row reduction (RREF)
A_extended = Matrix([[3, 2, -1, 1],
                     [2, -2, 4, -2],
                     [-1, 0.5, -1, 0]])

rref_form, pivot_cols = A_extended.rref()
print(f"RREF:\n{rref_form}")
print(f"Pivot columns: {pivot_cols}")

特殊函数和常数

数学常数

from sympy import pi, E, I, oo, zoo
from sympy import GoldenRatio, EulerGamma
from sympy import S  # Singleton for exact numbers

# Pi
print(f"π = {pi}")
print(f"π ≈ {pi.evalf()}")  # Numerical value

# Euler's number
print(f"e = {E}")
print(f"e ≈ {E.evalf()}")

# Imaginary unit
print(f"i = {I}")
print(f"i² = {I**2}")  # -1

# Infinity
print(f"∞ = {oo}")
print(f"Complex infinity = {zoo}")

# Golden ratio
print(f"φ = {GoldenRatio}")
print(f"φ ≈ {GoldenRatio.evalf()}")

# Euler-Mascheroni constant
print(f"γ = {EulerGamma}")
print(f"γ ≈ {EulerGamma.evalf()}")

# Exact rational numbers
from sympy import Rational
half = Rational(1, 2)
third = Rational(1, 3)
print(f"1/2 + 1/3 = {half + third}")  # 5/6

# Infinity arithmetic
print(f"1/∞ = {1/oo}")  # 0
print(f"∞ + 1 = {oo + 1}")  # oo

特殊函数

from sympy import symbols, factorial, binomial, fibonacci
from sympy import gamma, beta, zeta
from sympy import besselj, bessely

n, m, x = symbols('n m x')

# Factorial
print(f"5! = {factorial(5)}")  # 120
print(f"n! = {factorial(n)}")

# Binomial coefficient
print(f"C(5,2) = {binomial(5, 2)}")  # 10
print(f"C(n,m) = {binomial(n, m)}")

# Fibonacci
print(f"F(10) = {fibonacci(10)}")  # 55

# Gamma function
print(f"Γ(5) = {gamma(5)}")  # 24 = 4!
print(f"Γ(1/2) = {gamma(S(1)/2)}")  # sqrt(pi)

# Beta function
print(f"B(2,3) = {beta(2, 3)}")

# Riemann zeta function
print(f"ζ(2) = {zeta(2)}")  # pi**2/6
print(f"ζ(2) ≈ {zeta(2).evalf()}")

# Bessel functions
print(f"J₀(x) = {besselj(0, x)}")
print(f"J₁(x) = {besselj(1, x)}")

误差函数和超几何函数

from sympy import symbols, erf, erfc, erfi
from sympy import hyper, meijerg
from sympy import exp, sqrt, pi

x = symbols('x')

# Error function
print(f"erf(x) = {erf(x)}")
print(f"erf(1) ≈ {erf(1).evalf()}")

# Complementary error function
print(f"erfc(x) = {erfc(x)}")
print(f"erfc(x) = {1 - erf(x)}")  # Relationship

# Imaginary error function
print(f"erfi(x) = {erfi(x)}")

# Integral representation
from sympy import integrate
erf_integral = 2/sqrt(pi) * integrate(exp(-x**2), (x, 0, x))
print(f"erf via integral: {erf_integral}")

# Hypergeometric function
# 2F1(a, b; c; z)
result = hyper([1, 2], [3], x)
print(f"Hypergeometric: {result}")

代码生成

Lambdify（转换为 NumPy）

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

x = symbols('x')

# Symbolic expression
expr = sin(x) * exp(-x**2/2)

# Convert to numerical function
f = lambdify(x, expr, 'numpy')

# Evaluate on array
x_vals = np.linspace(-3, 3, 1000)
y_vals = f(x_vals)

# Plot
plt.figure(figsize=(10, 6))
plt.plot(x_vals, y_vals)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title(r'$f(x) = \sin(x) \cdot e^{-x^2/2}$')
plt.grid(True)
plt.savefig('lambdify_example.png', dpi=150)

print(f"Function evaluated at x=1: {f(1):.6f}")

# Multiple arguments
x, y = symbols('x y')
expr = x**2 + y**2
f = lambdify((x, y), expr, 'numpy')

X, Y = np.meshgrid(np.linspace(-2, 2, 50),
                    np.linspace(-2, 2, 50))
Z = f(X, Y)

print(f"Shape of result: {Z.shape}")

LaTeX 生成

from sympy import symbols, sin, cos, exp, integrate, Integral
from sympy import latex, Derivative, limit

x, y = symbols('x y')

# Convert expression to LaTeX
expr = sin(x)**2 + cos(x)**2
latex_str = latex(expr)
print(f"LaTeX: {latex_str}")

# Integral in LaTeX
expr = integrate(x**2, x)
latex_str = latex(expr)
print(f"∫x²dx in LaTeX: {latex_str}")

# Unevaluated integral
expr = Integral(sin(x) * exp(-x), (x, 0, float('inf')))
latex_str = latex(expr)
print(f"Integral notation: {latex_str}")

# Derivative in LaTeX
expr = Derivative(sin(x**2), x)
latex_str = latex(expr)
print(f"Derivative: {latex_str}")

# Complex expression
expr = limit((1 + 1/x)**x, x, float('inf'))
latex_str = latex(expr)
print(f"Limit: {latex_str}")

# Matrix in LaTeX
from sympy import Matrix
M = Matrix([[x, y], [y, x]])
latex_str = latex(M)
print(f"Matrix:\n{latex_str}")

代码生成（C、Fortran）

from sympy import symbols, sin, cos, exp
from sympy.utilities.codegen import codegen

x, y = symbols('x y')

# Expression to convert
expr = sin(x) * exp(-y)

# Generate C code
[(c_name, c_code), (h_name, c_header)] = codegen(
    ('func', expr), 'C', 'file', header=False
)

print("C code:")
print(c_code)

# Generate Fortran code
[(f_name, f_code)] = codegen(
    ('func', expr), 'F95', 'file', header=False
)

print("\nFortran code:")
print(f_code)

# More complex: function with multiple outputs
from sympy import symbols, Function, Derivative
from sympy.utilities.autowrap import autowrap

x = symbols('x')
expr1 = x**2
expr2 = x**3

# Can use autowrap for automatic compilation
# (requires compiler)
# f = autowrap(expr1, language='C', backend='Cython')

SymPy 转 NumPy 函数

from sympy import symbols, Matrix, lambdify
import numpy as np

x, y = symbols('x y')

# Matrix operations
M = Matrix([[x, y], [y, x]])

# Convert matrix to function
M_func = lambdify((x, y), M, 'numpy')

# Evaluate
result = M_func(1, 2)
print(f"Matrix at (1, 2):\n{result}")

# Matrix with operations
expr = M.det()  # Determinant
det_func = lambdify((x, y), expr, 'numpy')

print(f"Determinant at (1, 2): {det_func(1, 2)}")

# Vector operations
v = Matrix([x**2, y**2, x*y])
v_func = lambdify((x, y), v, 'numpy')

result = v_func(np.array([1, 2, 3]), np.array([4, 5, 6]))
print(f"Vector result:\n{result}")

物理与工程

经典力学

from sympy import symbols, Function, diff, solve, Eq
from sympy import sin, cos, sqrt, simplify

t = symbols('t', positive=True)
m, g, L = symbols('m g L', positive=True)

# Simple pendulum
theta = Function('theta')

# Equation of motion: theta'' + (g/L)*sin(theta) = 0
eq = Eq(diff(theta(t), t, 2) + (g/L)*sin(theta(t)), 0)

print(f"Pendulum equation: {eq}")

# Small angle approximation: sin(theta) ≈ theta
eq_linear = Eq(diff(theta(t), t, 2) + (g/L)*theta(t), 0)

# Solve linear equation
from sympy import dsolve
solution = dsolve(eq_linear, theta(t))
print(f"Small angle solution: {solution}")

# Harmonic oscillator
x = Function('x')
k = symbols('k', positive=True)

# mx'' + kx = 0
eq = Eq(m*diff(x(t), t, 2) + k*x(t), 0)
solution = dsolve(eq, x(t))
print(f"Harmonic oscillator: {solution}")

# Angular frequency
omega = sqrt(k/m)
print(f"Angular frequency: ω = {omega}")

量子力学

from sympy import symbols, Function, diff, exp, sqrt, integrate
from sympy import pi, I, oo, simplify
from sympy import Rational

x, t = symbols('x t', real=True)
m, hbar, omega = symbols('m hbar omega', positive=True)

# Quantum harmonic oscillator wavefunction (ground state)
psi_0 = (m*omega/(pi*hbar))**Rational(1,4) * exp(-m*omega*x**2/(2*hbar))

print(f"Ground state: ψ₀(x) = {psi_0}")

# Verify normalization
norm = integrate(psi_0**2, (x, -oo, oo))
norm_simplified = simplify(norm)
print(f"Normalization: ∫|ψ₀|²dx = {norm_simplified}")

# Expectation value of x
expect_x = integrate(x * psi_0**2, (x, -oo, oo))
print(f"⟨x⟩ = {simplify(expect_x)}")

# Expectation value of x²
expect_x2 = integrate(x**2 * psi_0**2, (x, -oo, oo))
expect_x2_simplified = simplify(expect_x2)
print(f"⟨x²⟩ = {expect_x2_simplified}")

# Uncertainty
delta_x = sqrt(expect_x2_simplified)
print(f"Δx = {delta_x}")

# Momentum operator expectation
# ⟨p⟩ = -iℏ∫ ψ* (dψ/dx) dx
dpsi_dx = diff(psi_0, x)
expect_p = integrate(-I*hbar*psi_0*dpsi_dx, (x, -oo, oo))
print(f"⟨p⟩ = {simplify(expect_p)}")

电磁学

from sympy import symbols, sin, cos, exp, sqrt, pi
from sympy import Function, diff, integrate
from sympy import I, oo

x, y, z, t = symbols('x y z t', real=True)
c, epsilon_0, mu_0 = symbols('c epsilon_0 mu_0', positive=True)

# Maxwell's equations in vacuum (symbolic)
E = Function('E')  # Electric field
B = Function('B')  # Magnetic field

# Gauss's law: ∇·E = 0 (in vacuum)
div_E = diff(E(x, y, z, t), x) + diff(E(x, y, z, t), y) + diff(E(x, y, z, t), z)

# Plane wave solution
k, omega = symbols('k omega', real=True)
E_0 = symbols('E_0', real=True)

# Electric field: E = E₀ exp(i(kx - ωt))
E_wave = E_0 * exp(I*(k*x - omega*t))
B_wave = E_wave / c  # Magnetic field

print(f"Electric field: E = {E_wave}")
print(f"Magnetic field: B = {B_wave}")

# Dispersion relation
dispersion = omega - c*k
print(f"Dispersion: ω = {c*k}")

# Poynting vector: S = (1/μ₀) E × B
# For plane wave in x-direction
S = E_0**2 / (mu_0 * c)
print(f"Time-averaged Poynting vector: ⟨S⟩ = {S}")

# Energy density
u = epsilon_0 * E_0**2 / 2
print(f"Energy density: u = {u}")

数论与离散数学

素数与因式分解

from sympy import primefactors, factorint, isprime, prime
from sympy import nextprime, prevprime, primerange
from sympy import divisors, totient

# Check if prime
n = 17
print(f"{n} is prime: {isprime(n)}")

# Prime factorization
n = 360
factors = factorint(n)
print(f"Prime factorization of {n}: {factors}")
# {2: 3, 3: 2, 5: 1} means 2³ · 3² · 5

# List of prime factors
primes = primefactors(n)
print(f"Prime factors: {primes}")

# n-th prime
print(f"10th prime: {prime(10)}")  # 29

# Next/previous prime
p = 100
print(f"Next prime after {p}: {nextprime(p)}")
print(f"Previous prime before {p}: {prevprime(p)}")

# Primes in range
primes_list = list(primerange(1, 50))
print(f"Primes up to 50: {primes_list}")

# Divisors
divs = divisors(n)
print(f"Divisors of {n}: {divs}")

# Euler's totient function
phi = totient(n)
print(f"φ({n}) = {phi}")

模运算

from sympy import symbols, Mod, gcd, lcm
from sympy import mod_inverse, crt, isprime

# Modular arithmetic
a, b, m = symbols('a b m', integer=True)

# Basic modular operations
print(f"17 mod 5 = {Mod(17, 5)}")  # 2

# GCD and LCM
print(f"gcd(48, 18) = {gcd(48, 18)}")  # 6
print(f"lcm(48, 18) = {lcm(48, 18)}")  # 144

# Modular inverse
# Find x such that (a*x) mod m = 1
a, m = 3, 11
inv = mod_inverse(a, m)
print(f"{a}⁻¹ mod {m} = {inv}")
# Verify
print(f"Verification: {Mod(a * inv, m)}")  # Should be 1

# Chinese Remainder Theorem
# Solve: x ≡ 2 (mod 3), x ≡ 3 (mod 5), x ≡ 2 (mod 7)
remainders = [2, 3, 2]
moduli = [3, 5, 7]
solution = crt(moduli, remainders)
print(f"CRT solution: x = {solution[0]} (mod {solution[1]})")

# Verify
x = solution[0]
for r, m in zip(remainders, moduli):
    print(f"x mod {m} = {Mod(x, m)} (should be {r})")

组合数学

from sympy import symbols, factorial, binomial, multinomial
from sympy import catalan, bell, fibonacci, lucas
from sympy import stirling, partition

n, k = symbols('n k', integer=True, positive=True)

# Factorials
print(f"10! = {factorial(10)}")

# Binomial coefficients
print(f"C(10, 3) = {binomial(10, 3)}")  # 120

# Multinomial coefficients
# (n choose k1, k2, ..., km) = n! / (k1! k2! ... km!)
print(f"Multinomial(5, [2, 2, 1]) = {multinomial(5, [2, 2, 1])}")

# Catalan numbers
for i in range(6):
    print(f"C_{i} = {catalan(i)}")

# Bell numbers (number of partitions of a set)
for i in range(6):
    print(f"B_{i} = {bell(i)}")

# Fibonacci and Lucas numbers
print(f"Fibonacci: {[fibonacci(i) for i in range(10)]}")
print(f"Lucas: {[lucas(i) for i in range(10)]}")

# Stirling numbers (first kind)
print(f"s(5, 2) = {stirling(5, 2, kind=1)}")

# Stirling numbers (second kind)
print(f"S(5, 2) = {stirling(5, 2, kind=2)}")

# Integer partitions
from sympy.utilities.iterables import partitions
p = list(partitions(5))
print(f"Partitions of 5: {p}")

实用工作流

符号回归分析

from sympy import symbols, exp, log, sin, cos
from sympy import lambdify, simplify
import numpy as np
from scipy.optimize import curve_fit

# Generate synthetic data
x_data = np.linspace(0, 10, 50)
y_true = 2.5 * np.exp(-0.5 * x_data) * np.sin(x_data)
y_data = y_true + 0.1 * np.random.randn(len(x_data))

# Guess functional form
x = symbols('x')
A, alpha, omega = symbols('A alpha omega', real=True)

# Model: y = A * exp(-alpha*x) * sin(omega*x)
model = A * exp(-alpha*x) * sin(omega*x)

print(f"Model: y = {model}")

# Convert to numerical function
f_model = lambdify((x, A, alpha, omega), model, 'numpy')

# Fit using scipy
def fit_func(x_vals, A_val, alpha_val, omega_val):
    return f_model(x_vals, A_val, alpha_val, omega_val)

params, covariance = curve_fit(fit_func, x_data, y_data, p0=[2, 0.5, 1])

A_fit, alpha_fit, omega_fit = params
print(f"\nFitted parameters:")
print(f"A = {A_fit:.4f} (true: 2.5)")
print(f"α = {alpha_fit:.4f} (true: 0.5)")
print(f"ω = {omega_fit:.4f} (true: 1.0)")

# Create fitted expression
fitted_expr = model.subs([(A, A_fit), (alpha, alpha_fit), (omega, omega_fit)])
print(f"\nFitted expression: y = {fitted_expr}")

# Plot
import matplotlib.pyplot as plt

y_fit = f_model(x_data, A_fit, alpha_fit, omega_fit)

plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, alpha=0.5, label='Data')
plt.plot(x_data, y_true, 'g--', label='True')
plt.plot(x_data, y_fit, 'r-', label='Fit')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.savefig('symbolic_fit.png', dpi=150)

带约束的优化问题

from sympy import symbols, diff, solve, simplify
from sympy import Matrix, hessian

# Optimization: minimize f(x,y) = x² + y²
# Subject to: g(x,y) = x + y - 1 = 0

x, y, lam = symbols('x y lambda', real=True)

# Objective function
f = x**2 + y**2

# Constraint
g = x + y - 1

# Lagrangian: L = f - λg
L = f - lam * g

# KKT conditions: ∇L = 0
dL_dx = diff(L, x)
dL_dy = diff(L, y)
dL_dlam = diff(L, lam)

print("KKT conditions:")
print(f"∂L/∂x = {dL_dx} = 0")
print(f"∂L/∂y = {dL_dy} = 0")
print(f"∂L/∂λ = {dL_dlam} = 0")

# Solve system
solution = solve([dL_dx, dL_dy, dL_dlam], [x, y, lam])
print(f"\nSolution: {solution}")

x_opt = solution[x]
y_opt = solution[y]

print(f"\nOptimal point: ({x_opt}, {y_opt})")
print(f"Objective value: f({x_opt}, {y_opt}) = {f.subs([(x, x_opt), (y, y_opt)])}")

# Verify constraint
g_check = g.subs([(x, x_opt), (y, y_opt)])
print(f"Constraint check: g = {simplify(g_check)}")

# Second-order conditions (Hessian)
H = hessian(f, (x, y))
print(f"\nHessian of f:\n{H}")

# Evaluate at optimum
H_opt = H.subs([(x, x_opt), (y, y_opt)])
eigenvals = H_opt.eigenvals()
print(f"Eigenvalues at optimum: {eigenvals}")
if all(val > 0 for val in eigenvals.keys()):
    print("✓ Confirmed: Local minimum")

符号微分方程求解器

from sympy import symbols, Function, dsolve, Eq, diff
from sympy import sin, cos, exp, sqrt
from sympy import init_printing

init_printing()

x = symbols('x')
y = Function('y')

# Collection of ODEs to solve
odes = [
    # Linear first-order
    (Eq(diff(y(x), x) + y(x), exp(x)), "Linear first-order"),
    
    # Separable
    (Eq(diff(y(x), x), y(x) * (1 - y(x))), "Logistic equation"),
    
    # Second-order linear
    (Eq(diff(y(x), x, 2) - 4*y(x), 0), "Simple harmonic"),
    
    # Damped oscillator
    (Eq(diff(y(x), x, 2) + 2*diff(y(x), x) + 5*y(x), 0), "Damped oscillator"),
]

print("SOLVING DIFFERENTIAL EQUATIONS")
print("=" * 60)

for eq, description in odes:
    print(f"\n{description}:")
    print(f"Equation: {eq}")
    
    # Solve
    solution = dsolve(eq, y(x))
    print(f"Solution: {solution}")
    
    # Verify solution
    y_sol = solution.rhs
    eq_check = eq.lhs.subs(y(x), y_sol)
    eq_check_simplified = simplify(eq_check)
    
    if eq_check_simplified == eq.rhs:
        print("✓ Solution verified")
    else:
        print(f"Verification: {eq_check_simplified} = {eq.rhs}")

# Solve with initial conditions
print("\n" + "=" * 60)
print("WITH INITIAL CONDITIONS")
print("=" * 60)

# y' + y = 0, y(0) = 1
eq = Eq(diff(y(x), x) + y(x), 0)
ics = {y(0): 1}

solution = dsolve(eq, y(x), ics=ics)
print(f"Equation: {eq}, y(0) = 1")
print(f"Solution: {solution}")

泰勒级数近似分析

from sympy import symbols, sin, cos, exp, log, sqrt
from sympy import series, diff, factorial
from sympy import lambdify, Abs
import numpy as np
import matplotlib.pyplot as plt

x = symbols('x')

# Function to approximate
f = exp(sin(x))

# Calculate Taylor series of different orders
orders = [1, 3, 5, 7, 9]
series_dict = {}

print(f"Taylor series for f(x) = {f} around x=0:")
print("=" * 60)

for n in orders:
    s = series(f, x, 0, n)
    s_poly = s.removeO()
    series_dict[n] = s_poly
    print(f"Order {n-1}: {s}")

# Convert to numerical functions
f_num = lambdify(x, f, 'numpy')
series_num = {n: lambdify(x, s, 'numpy') for n, s in series_dict.items()}

# Evaluate
x_vals = np.linspace(-2, 2, 200)
y_true = f_num(x_vals)

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Series approximations
ax1.plot(x_vals, y_true, 'k-', linewidth=2, label='True function')
for n in orders:
    y_approx = series_num[n](x_vals)
    ax1.plot(x_vals, y_approx, '--', label=f'Order {n-1}', alpha=0.7)

ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_title('Taylor Series Approximations')
ax1.legend()
ax1.grid(True)
ax1.set_ylim(-2, 5)

# Error analysis
for n in orders:
    y_approx = series_num[n](x_vals)
    error = np.abs(y_true - y_approx)
    ax2.semilogy(x_vals, error, label=f'Order {n-1}')

ax2.set_xlabel('x')
ax2.set_ylabel('Absolute Error')
ax2.set_title('Approximation Error')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.savefig('taylor_series_analysis.png', dpi=150)

# Print error at specific points
print("\n" + "=" * 60)
print("ERROR AT x = 1:")
print("=" * 60)
x_test = 1.0
y_test = f_num(x_test)

for n in orders:
    y_approx = series_num[n](x_test)
    error = abs(y_test - y_approx)
    rel_error = error / abs(y_test) * 100
    print(f"Order {n-1}: Error = {error:.6e} ({rel_error:.2f}%)")

符号积分技巧

from sympy import symbols, integrate, sin, cos, exp, log, sqrt
from sympy import apart, trigsimp, simplify
from sympy import pi, oo, I

x = symbols('x')

# Collection of integrals demonstrating different techniques
integrals = [
    # Substitution (automatic)
    (x * exp(x**2), "u-substitution (automatic)"),
    
    # Integration by parts (automatic)
    (x * sin(x), "Integration by parts"),
    
    # Trigonometric
    (sin(x)**2, "Trigonometric identity"),
    
    # Rational functions (partial fractions)
    (1/(x**2 + 3*x + 2), "Partial fractions"),
    
    # Logarithmic
    (log(x)/x, "Logarithmic"),
    
    # Radical
    (1/sqrt(1 - x**2), "Inverse trig"),
]

print("INTEGRATION TECHNIQUES")
print("=" * 60)

for integrand, technique in integrals:
    print(f"\n{technique}:")
    print(f"∫ {integrand} dx")
    
    # Compute integral
    result = integrate(integrand, x)
    print(f"= {result}")
    
    # Verify by differentiation
    check = diff(result, x)
    check_simplified = simplify(check)
    
    if simplify(check_simplified - integrand) == 0:
        print("✓ Verified by differentiation")
    else:
        print(f"Check: d/dx({result}) = {check_simplified}")

# Definite integrals
print("\n" + "=" * 60)
print("DEFINITE INTEGRALS")
print("=" * 60)

definite_integrals = [
    (x**2, 0, 1, "∫₀¹ x² dx"),
    (sin(x), 0, pi, "∫₀^π sin(x) dx"),
    (exp(-x**2), -oo, oo, "∫₋∞^∞ e^(-x²) dx"),
    (1/(1 + x**2), -oo, oo, "∫₋∞^∞ 1/(1+x²) dx"),
]

for integrand, a, b, description in definite_integrals:
    result = integrate(integrand, (x, a, b))
    print(f"{description} = {result}")
    if result.is_number:
        print(f"  ≈ {result.evalf()}")

这份全面的 SymPy 指南涵盖了 60 多个示例，涉及符号数学、微积分、代数等内容！