Python的符号计算：Sympy

samonyu
201 阅读8分钟

该文章为《Numerical Python 第三版》第三章符号计算的内容的整理。

符号计算

符号计算是一种与基于数值数组的计算完全不同的范式。在符号计算软件中，也称为计算机代数系统（CAS），数学对象和表达式的表示被分析地操作和变换。符号计算主要是使用计算机自动化手工纸笔可以完成的分析计算。然而，通过使用计算机代数系统自动化书写的记录和数学表达式的操作，可以使分析计算远远超越手工所能达到的范围。

分析与符号计算是科学和技术计算领域的重要组成部分。即使对于只能通过数值方法解决的问题（由于许多实际问题中分析方法不可行），在不得不采用数值技术之前尽可能利用分析方法所能做到的事情可以有很大的不同。例如：

  • 减少数值问题的复杂性或规模：通过分析方法简化问题，可以降低最终需要解决的数值问题的复杂度，从而提高计算效率。
  • 提供数值计算的初始条件或边界条件：符号计算可以为数值计算提供精确的初始条件或边界条件，从而提高数值计算的精度。

SymPy：Python 中的符号计算库

在科学计算的 Python 环境中，主要的符号计算库是 SymPy（Symbolic Python）。SymPy 是用 Python 编写的，并提供了广泛用于分析和符号问题的工具。

SymPy 的主要功能

  • 符号表达式的创建和操作：SymPy 可以创建符号变量、函数、表达式，并进行各种代数运算，例如展开、因式分解、化简等。
  • 微积分：SymPy 可以进行符号微分、积分、极限计算等。
  • 方程求解：SymPy 可以求解代数方程、微分方程等。
  • 线性代数：SymPy 可以进行矩阵运算、行列式计算、特征值求解等。
  • 绘图：SymPy 可以绘制符号表达式的图形。

Sympy 导入

# 导入
# 通常的导入方法
from sympy import *
# 避免冲突的导入方法
import sympy

# 配置打印系统
sympy.init_printing()

# 导入常用符号（不是数值）
from sympy import I, pi, oo

符号

# 创建一个符号（字符串）
x = sympy.Symbol("x")

# 创建一个实数符号
y = sympy.Symbol("y", real=True)

# 开方
# 不可以简化
x = sympy.Symbol("x")
sympy.sqrt(x ** 2)
# 可以简化（进一步计算）
y = sympy.Symbol("y", positive=True)
sympy.sqrt(y ** 2)

# 创建多个符号
a, b, c = sympy.symbols("a, b, c", negative=True)

符号的属性

关键词参数属性描述
​real​, imaginary​​is_real​, is_imaginary​符号代表实数或虚数
​positive​, negative​​is_positive​, is_negative​符号代表整数或复数
​Integer​​is_integer​符号代表一个整数
​odd​, even​​is_odd​, is_even​符号代表奇数或偶数
​Prime​​is_prime​符号代表一个质数（也是整数）
​finite​, infinite​​is_finite​, is_infinite​符号代表有限或无限量

数字

# 创建数字
i = sympy.Integer(19)
f = sympy.Float(2.3)
i, f = sympy.sympify(19), sympy.sympify(2.3)

# 转换为Python内置类型
int(i)
float(f)

# 处理特别大的数
# 100的阶乘
sympy.factorial(100)

# 精确的浮点数，使用字符串初始化
sympy.Float('0.3', 25)

# 有理数（分数）
sympy.Rational(11, 13)

特殊的常数

数学符号SymPy 符号描述
π​sympy.pi​圆的周长与直径的比值
e​sympy.E​自然对数的底数，​
γ​sympy.EulerGamma​欧拉常数
i​sympy.I​虚数单位，​
∞​sympy.oo​无穷大

函数

# 声明函数
x, y, z = sympy.symbols("x, y, z")
f = sympy.Function("f")
g = sympy.Function("g")(x, y, z)

# lambda 函数
h = sympy.Lambda(x, x**2)

表达式

x = sympy.Symbol("x")
# 表达式
expr = 1 + 2 * x**2 + 3 * x**3

​image转存失败，建议直接上传图片文件​

对表达式进行操作

化简

SymPy 中的表达式应被视为不可变对象（不能或不应该就地更改）。这些函数通常不会改变传递给它们的表达式，而是创建一个新的表达式来表示修改后的形式。

函数描述
​sympy.simplify​尝试各种方法和途径以获得给定表达式的简化形式。
​sympy.trigsimp​尝试使用三角恒等式简化表达式。
​sympy.powsimp​尝试使用幂的法则简化表达式。
​sympy.compsimp​简化组合表达式。
​sympy.ratsimp​通过将表达式写在一个共同的分母上来简化表达式。
expr = 2 * (x**2 - x) - x * (x + 1)
# OUT: 2x2 – x(x+1)–2x

# 化简表达式
sympy.simplify(expr)
# OUT: x(x–3)

# 化简表达式
expr.simplify()
# OUT: x(x–3)

# 原表达式不改变
expr
# OUT: 2x2 – x(x+1)–2x

展开

函数描述参数
​sympy.expand​对表达式进行展开展开乘积：mul=True​ 展开三角函数：trig=True​ 展开对数：log=True​ 展开实部与虚部：complex=True​ 展开底数：power_base=True​ 展开指数：power_exp=True​

因式分解、提取

# 对x进行因式分解
sympy.factor(x**2 - 1)
# OUT: (x – 1)(x + 1)

# 因式分解
sympy.factor(x * sympy.cos(y) + sympy.sin(z) * x)
# OUT: x(sin(z) + cos(y))

# 对数提取
sympy.logcombine(sympy.log(a) - sympy.log(b))

# 对x、y提取（结果不同）
expr = x + y + x * y * z
expr.collect(x)
expr.collect(y)

分数的处理

# 分数的拆分
sympy.apart(1/(x**2 + 3*x + 2), x)

# 分数的合并
sympy.together(1 / (y * x + y) + 1 / (1+x))

# 分数的消去共享因子
sympy.cancel(y / (y * x + y))

替换

# 将x替换为y
(x + y).subs(x, y)
sympy.sin(x * sympy.exp(x)).subs(x, y)

# 多个替换
sympy.sin(x * z).subs({z: sympy.exp(y), x: y, sympy.sin: sympy.cos})

# 符号替换成值（代入值）
expr = x * y + z**2 *x
values = {x: 1.25, y: 0.4, z: 3.2} 
expr.subs(values)

数值评估

一个 SymPy 表达式可以通过 sympy.N​ 函数或 SymPy 表达式实例的 evalf​ 方法来评估。

无论是 sympy.N​ 函数还是 evalf​ 方法，都接受一个可选参数，用于指定表达式计算时的有效位数。

当我们需要对一系列输入值进行数值评估时，可以通过遍历这些值并依次调用 evalf​ 方法来实现。

​sympy.lambdify​ 函数接受一组自由符号和一个表达式作为参数，并生成一个高效计算表达式数值的函数。

sympy.N(1 + pi)

# 有效位数
sympy.N(pi, 50)
(x + 1/pi).evalf(10)

# 对一系列输入值进行评估
expr = sympy.sin(pi * x * sympy.exp(x))
[expr.subs(x, xx).evalf(3) for xx in range(0, 10)]

# 更快的评估
expr_func = sympy.lambdify(x, expr)
expr_func(1.0)
[expr_func(i) for i in range (0,10)]

# 与numpy结合
expr_func = sympy.lambdify(x, expr, 'numpy')
import numpy as np
xvalues = np.arange(0, 10)
expr_func(xvalues)

微积分

微分

# 一元微分（一阶、二阶、三阶）
f = sympy.Function('f')(x)
sympy.diff(f, x)		# equivalent to f.diff(x)
sympy.diff(f, x, x)
sympy.diff(f, x, 3)		# equivalent to sympy.diff(f, x, x, x)

# 多元微分
g = sympy.Function('g')(x, y)
g.diff(x, y)			# equivalent to sympy.diff(g, x, y)
g.diff(x, 3, y, 2)    	# equivalent to sympy.diff(g, x, x, x, y, y)

# 显式的求值
.doit()

积分

a, b, x, y = sympy.symbols("a, b, x, y")
f = sympy.Function("f")(x)
# 对x进行积分
sympy.integrate(f)
# 对x从a到b进行积分
sympy.integrate(f, (x, a, b))

# 多变量的积分，显式指定积分变量
expr = sympy.sin(x*sympy.exp(y))
sympy.integrate(expr, x)

# 多重积分
expr = (x + y)**2
sympy.integrate(expr, (x, 0, 1), (y, 0, 1))

级数展开

在SymPy中，可以通过函数sympy.series​或SymPy表达式实例中的series​方法计算函数或表达式的级数展开。

​sympy.series​的第一个参数是要展开的函数或表达式，随后是一个用以计算展开的符号（对于单变量的表达式和函数，该参数可以省略）。

此外，还可以指定级数展开的具体点（使用x0​关键字参数，默认x0=0​），指定展开的阶数（使用n​关键字参数，默认n=6​），以及指定从x0​的上下方向中何处开始计算级数（使用dir​关键字参数，默认dir='+'​）。

x, y = sympy.symbols("x, y")
f = sympy.Function("f")(x)
# 在0点，展开6阶，正向
sympy.series(f, x)

# 在x=x0点，展开2阶
x0 = sympy.Symbol("{x_0}")
f.series(x, x0, n=2)
# 去除小o
f.series(x, x0, n=2).removeO()

极限

在 SymPy 中，可以使用 sympy.limit​ 函数来计算极限。该函数接受一个表达式、依赖的符号以及符号在极限中的趋近值。

sympy.limit(sympy.sin(x) / x, x, 0)

求和与乘积

n = sympy.symbols("n", integer=True)

x = sympy.Sum(1/(n**2), (n, 1, oo))
x.doit()

x = sympy.Product(n, (n, 1, 7))
x.doit()

方程

x = sympy.Symbol("x")
# 解方程=0，求x
sympy.solve(x**2 + 2*x - 3)
# OUT: [-3, 1]

# 一元二次方程通解
a, b, c = sympy.symbols("a, b, c")
sympy.solve(a * x**2 + b * x + c, x)

# 解方程组
eq1 = x + 2 * y – 1
eq2 = x - y + 1
# 每个解以字典的方式返回
sympy.solve([eq1, eq2], [x, y], dict=True)

线性代数

a, b, c, d = sympy.symbols("a, b, c, d")
M = sympy.Matrix([[a, b], [c, d]])

sympy.Matrix(3, 4, lambda m, n: 10 * m + n)
函数/方法描述
​M.transpose()​/M.T​计算矩阵的转置。
​M.adjoint()​/M.H​计算矩阵的伴随矩阵。
​M.trace()​计算矩阵的迹（对角线元素之和）。
​M.det()​计算矩阵的行列式。
​M.inv()​计算矩阵的逆矩阵。
​M.LUdecomposition()​计算矩阵的 LU 分解。
​M.LUsolve(b)​使用 LU 分解求解线性方程组 Mx = b​，其中 x​ 为未知向量。
​M.QRdecomposition()​计算矩阵的 QR 分解。
​M.QRsolve(b)​使用 QR 分解求解线性方程组 Mx = b​，其中 x​ 为未知向量。
​M.diagonalize()​对角化矩阵 M​，使其可以表示为 D = P⁻¹MP​ 的形式，其中 D​ 为对角矩阵。
​M.norm()​计算矩阵的范数。
​M.nullspace()​计算矩阵的零空间的一组基向量。
​M.rank()​计算矩阵的秩。
​M.singular_values()​计算矩阵的奇异值。
​M.solve(b)​求解线性方程组 Mx = b​。

参考文献

  1. R. Johansson, Numerical Python, doi.org/10.1007/979…
  2. www.sympy.org/en/index.ht…
avatar
文章
阅读
粉丝