SymPy-1-13-中文文档-二-
SymPy 1.13 中文文档（二）
指南
指南是关于如何执行特定任务的逐步说明。
对于其他 SymPy 主题的深入详尽探索，请参阅 解释 和 API 参考 部分。
-
假设
-
符号和模糊布尔值
-
编写自定义函数
-
物理学
-
解方程
-
SymPy 标志
-
引用 SymPy
假设
本页概述了 SymPy 中的核心假设系统。它解释了核心假设系统是什么，假设系统如何使用以及不同的假设谓词的含义。
注意
本页描述了核心假设系统，通常也称为“旧假设”系统。还有一个“新假设”系统，描述在其他地方。请注意，这里描述的系统实际上是 SymPy 中广泛使用的系统。“新假设”系统目前实际上在 SymPy 中没有得到使用，而“旧假设”系统也不会被移除。在撰写本文时（SymPy 1.7），仍建议用户使用旧的假设系统。
首先，我们考虑对具体整数（如 (2) 或 (-2)）进行平方根操作时会发生什么：
>>> from sympy import sqrt
>>> sqrt(2**2)
2
>>> sqrt((-2)**2)
2
>>> x = 2
>>> sqrt(x**2)
2
>>> sqrt(x**2) == x
True
>>> y = -2
>>> sqrt(y**2) == y
False
>>> sqrt(y**2) == -y
True
这些示例展示了对于正数 (x)，我们有 (\sqrt{x²} = x)，而对于负数，我们有 (\sqrt{x²} = -x)。这似乎是显而易见的，但在处理符号而非显式数字时，情况可能更令人惊讶。例如
>>> from sympy import Symbol, simplify
>>> x = Symbol('x')
>>> sqrt(x**2)
sqrt(x**2)
看起来应该简化为 x，但即使使用 simplify() 也不会：
>>> simplify(sqrt(x**2))
sqrt(x**2)
这是因为如果简化对于x的每一个可能值都不合适，SymPy 将拒绝简化这个表达式。默认情况下，符号 x 被认为仅代表大致像一个任意复数的东西，而这里的显而易见的简化仅对正实数有效。因为 x 不被认为是正数甚至是实数，所以这个表达式的简化是不可能的。
在创建符号时，我们可以告诉 SymPy 符号表示正实数，然后简化将自动发生。
>>> y = Symbol('y', positive=True)
>>> sqrt(y**2)
y
这就是 SymPy 中所谓的“假设”。如果用 positive=True 创建符号 y，那么 SymPy 将假设它表示一个正实数，而不是任意复数或可能是无限大的数。这种假设可以使表达式简化或可能允许其他操作。在创建符号时，尽可能明确符号的假设通常是一个好主意。
（旧）假设系统
假设系统有两个方面。一方面是我们可以在创建符号时声明符号的假设。另一方面是我们可以使用相应的 is_* 属性查询任何表达式的假设。例如：
>>> x = Symbol('x', positive=True)
>>> x.is_positive
True
我们可以查询任何表达式的假设，而不仅仅是一个符号：
>>> x = Symbol('x', positive=True)
>>> expr = 1 + x**2
>>> expr
x**2 + 1
>>> expr.is_positive
True
>>> expr.is_negative
False
假设查询中给出的值使用三值“模糊”逻辑。任何查询可以返回 True、False 或 None，其中 None 应被解释为结果是未知的。
>>> x = Symbol('x')
>>> y = Symbol('y', positive=True)
>>> z = Symbol('z', negative=True)
>>> print(x.is_positive)
None
>>> print(y.is_positive)
True
>>> print(z.is_positive)
False
注意
在上述示例中，我们需要使用 print，因为默认情况下 Python 解释器不显示特殊值 None。
有几个原因会导致假设查询返回 None。例如，查询可能是不可知的，就像上述的 x 情况一样。由于 x 没有声明任何假设，它大致代表一个任意复数。一个任意复数可能是正实数，但也可能不是。没有更多信息，无法解决查询 x.is_positive。
另一个假设查询返回 None 的原因是在许多情况下，例如确定表达式是否为正是不可判定的问题。这意味着通常情况下不存在回答该查询的算法。对于某些情况来说，可能存在算法或至少可以进行简单检查，但尚未在 SymPy 中实现，尽管可以添加。
假设查询可能返回 None 的最后一个原因是，假设系统并未尝试非常努力地回答复杂的查询。该系统旨在快速运行，并使用简单的启发式方法在常见情况下得出True或False的答案。例如，任何正项的和都是正的，因此：
>>> from sympy import symbols
>>> x, y = symbols('x, y', positive=True)
>>> expr = x + y
>>> expr
x + y
>>> expr.is_positive
True
最后一个示例特别简单，因此假设系统能够给出明确的答案。如果该和涉及正负项的混合，那么这将是一个更难的查询：
>>> x = Symbol('x', real=True)
>>> expr = 1 + (x - 2)**2
>>> expr
(x - 2)**2 + 1
>>> expr.is_positive
True
>>> expr2 = expr.expand()
>>> expr2
x**2 - 4*x + 5
>>> print(expr2.is_positive)
None
理想情况下，最后一个示例应该返回 True 而不是 None，因为该表达式对于任何实数值的 x 都是始终为正（且已假设 x 为实数）。尽管如此，假设系统的设计意在高效运行：预计许多更复杂的查询将无法完全解决。这是因为假设查询主要在 SymPy 内部作为低级别计算的一部分使用。使系统更全面可能会减慢 SymPy 的运行速度。
请注意，在模糊逻辑中，给出不确定结果 None 绝不是矛盾的。如果在解决查询时能推断出明确的 True 或 False 结果，则比返回 None 更好。然而，None 的结果并不是错误。任何使用假设系统的代码都需要准备处理任何查询的三种情况，并不应假设总会给出明确的答案。
假设系统不仅适用于符号或复杂表达式。它也可以用于普通的 SymPy 整数和其他对象。假设谓词适用于 Basic 的任何实例，这是 SymPy 对象大多数类的超类。普通的 Python int 不是 Basic 的实例，不能用于查询假设谓词。我们可以用 sympify() 或 S (SingletonRegistry) 将常规 Python 对象“符号化”为 SymPy 对象，然后可以使用假设系统：
>>> from sympy import S
>>> x = 2
>>> x.is_positive
Traceback (most recent call last):
...
AttributeError: 'int' object has no attribute 'is_positive'
>>> x = S(2)
>>> type(x)
<class 'sympy.core.numbers.Integer'>
>>> x.is_positive
True
注意：具有不同假设的符号
在 SymPy 中，可以声明两个具有不同名称的符号，它们将隐式地被视为 结构相等：
>>> x1 = Symbol('x')
>>> x2 = Symbol('x')
>>> x1
x
>>> x2
x
>>> x1 == x2
True
但是，如果符号具有不同的假设，则它们将被视为代表不同的符号：
>>> x1 = Symbol('x', positive=True)
>>> x2 = Symbol('x')
>>> x1
x
>>> x2
x
>>> x1 == x2
False
简化表达式的一种方法是使用 posify() 函数，它将表达式中的所有符号替换为具有假设 positive=True 的符号（除非这与符号的任何现有假设相矛盾）：
>>> from sympy import posify, exp
>>> x = Symbol('x')
>>> expr = exp(sqrt(x**2))
>>> expr
exp(sqrt(x**2))
>>> posify(expr)
(exp(_x), {_x: x})
>>> expr2, rep = posify(expr)
>>> expr2
exp(_x)
posify() 函数返回表达式，其中所有符号都被替换（这可能导致简化），还返回一个字典，将新符号映射到旧符号，可用于 subs()。这很有用，否则具有 positive=True 假设的新表达式将不等同于旧表达式：
>>> expr2
exp(_x)
>>> expr2 == exp(x)
False
>>> expr2.subs(rep)
exp(x)
>>> expr2.subs(rep) == exp(x)
True
对字符串输入应用假设
我们已经看到如何在 Symbol 或 symbols() 明确设置假设。一个自然的问题是，在什么其他情况下可以为对象分配假设？
用户通常将字符串作为 SymPy 函数的输入（尽管 SymPy 开发人员普遍认为这应该被弃用），例如：
>>> from sympy import solve
>>> solve('x**2 - 1')
[-1, 1]
当显式创建符号时，可以分配假设，这会影响 solve() 的行为：
>>> x = Symbol('x', positive=True)
>>> solve(x**2 - 1)
[1]
当使用字符串输入 SymPy 将创建表达式并隐式创建所有符号，因此问题是如何指定假设？答案是，与其依赖隐式字符串转换，不如显式使用parse_expr()函数，然后可以为符号提供假设，例如：
>>> from sympy import parse_expr
>>> parse_expr('x**2 - 1')
x**2 - 1
>>> eq = parse_expr('x**2 - 1', {'x':Symbol('x', positive=True)})
>>> solve(eq)
[1]
注意
solve()函数作为高级 API 是不寻常的，因为它实际上会检查输入符号（未知数）的假设，并使用它来定制其输出。假设系统否则影响低级评估，但高级 API 不一定显式处理。
谓词
有许多不同的谓词可以为符号假设或可以在表达式中查询。在创建符号时可以组合多个谓词。谓词在逻辑上使用and组合，因此如果一个符号被声明为positive=True并且也被声明为integer=True，那么它既是正数又是整数：
>>> x = Symbol('x', positive=True, integer=True)
>>> x.is_positive
True
>>> x.is_integer
True
可以使用assumptions0属性访问符号的完整已知谓词集：
>>> x.assumptions0
{'algebraic': True,
'commutative': True,
'complex': True,
'extended_negative': False,
'extended_nonnegative': True,
'extended_nonpositive': False,
'extended_nonzero': True,
'extended_positive': True,
'extended_real': True,
'finite': True,
'hermitian': True,
'imaginary': False,
'infinite': False,
'integer': True,
'irrational': False,
'negative': False,
'noninteger': False,
'nonnegative': True,
'nonpositive': False,
'nonzero': True,
'positive': True,
'rational': True,
'real': True,
'transcendental': False,
'zero': False}
我们可以看到列出的谓词比用于创建x的两个谓词多得多。这是因为假设系统可以从其他谓词的组合中推断出一些谓词。例如，如果一个符号声明为positive=True，那么可以推断它应该negative=False，因为正数永远不会是负数。类似地，如果一个符号被创建为integer=True，那么可以推断它应该rational=True，因为每个整数都是有理数。
下面提供了可能的谓词及其定义的完整表格。
假设（旧）假设的谓词
| 谓词 | 定义 | 含义 |
|---|---|---|
commutative |
可交换的表达式。 commutative表达式在乘法下与所有其他表达式交换。如果表达式a具有commutative=True，则对于任何其他表达式b都有a * b == b * a（即使b不是commutative）。与所有其他假设谓词不同，commutative必须始终是True或False，不能是None。而且与所有其他谓词不同，例如在Symbol('x')中，默认为True。[commutative] |
|
infinite |
无限表达式，如oo、-oo或zoo。[infinite] |
== !finite |
finite |
有限表达式。任何非infinite的表达式都被视为finite。[infinite] |
== !infinite |
hermitian |
Hermitian 操作符的域中的元素。[antihermitian] | |
antihermitian |
一个反厄米算子域的元素。[antihermitian] | |
complex |
复数，(z\in\mathbb{C})。任何形式为 (x + iy)，其中 (x) 和 (y) 是real且 (i = \sqrt{-1})。所有的complex数都是finite。包括所有的real数。 [complex] |
-> commutative``-> finite |
algebraic |
代数数，(z\in\overline{\mathbb{Q}})。任何是非零有理系数多项式 (p(z)\in\mathbb{Q}[z]) 的根的数。所有的algebraic数都是complex。一个algebraic数可能是real，也可能不是。包括所有的rational数。 [algebraic] |
-> complex |
transcendental |
复数，不属于代数数，(z\in\mathbb{C}-\overline{\mathbb{Q}})。所有的transcendental数都是complex。一个transcendental数可能是real，也可能不是，但绝不可能是rational。 [transcendental] |
== (complex & !algebraic) |
extended_real |
扩展实数线上的元素，(x\in\overline{\mathbb{R}}) 其中 (\overline{\mathbb{R}}=\mathbb{R}\cup{-\infty,+\infty})。一个extended_real数要么是real，要么是 (\pm\infty)。只有表达式为extended_real时才定义了关系运算符 <、<=、>= 和 >。 [extended_real] |
-> commutative |
real |
实数，(x\in\mathbb{R})。所有的real数都是finite和complex（实数集是复数集的子集）。包括所有的rational数。一个real数要么是negative、zero或positive。 [real] |
-> complex``== (extended_real & finite)``== (negative | zero | positive)``-> hermitian |
imaginary |
虚数，(z\in\mathbb{I}-{0})。形式为 (z=yi)，其中 (y) 是real，(y\ne 0) 且 (i=\sqrt{-1})。所有的imaginary数都是complex且不是real。特别地，SymPy 中不认为zero是imaginary。 [imaginary] |
-> complex``-> antihermitian``-> !extended_real |
rational |
有理数，(q\in\mathbb{Q})。任何形式为 (\frac{a}{b})，其中 (a) 和 (b) 是整数且 (b \ne 0)。所有的rational数都是real和algebraic。包括所有的integer数。 [rational] |
-> real``-> algebraic |
irrational |
不是有理数的实数，(x\in\mathbb{R}-\mathbb{Q})。[irrational] | == (real & !rational) |
integer |
整数，(a\in\mathbb{Z})。所有的整数都是rational。包括zero和所有的prime、composite、even和odd数。 [integer] |
-> rational |
noninteger |
不是整数的扩展实数，(x\in\overline{\mathbb{R}}-\mathbb{Z})。 | == (extended_real & !integer) |
even |
偶数，(e\in{2k: k\in\mathbb{Z}})。所有的 even 数都是整数。包括 zero。[奇偶性] |
-> 整数``-> !奇数 |
odd |
奇数，(o\in{2k + 1: k\in\mathbb{Z}})。所有的 odd 数都是整数。[奇偶性] |
-> 整数``-> !偶数 |
prime |
质数，(p\in\mathbb{P})。所有的 prime 数都是 positive 和 integer。[质数] |
-> 整数``-> 正数 |
composite |
复合数，(c\in\mathbb{N}-(\mathbb{P}\cup{1}))。是两个或更多素数的乘积的正整数。composite 数总是一个 positive integer，且不是 prime。[composite] |
-> (整数 & 正数 & !素数)``!composite -> (!正数 | !偶数 | 素数) |
zero |
数字 (0)。表达式中 zero=True 表示数字 0 是一个 integer。[零] |
-> 偶数 & 有限``== (扩展非负数 & 扩展非正数)``== (非负数 & 非正数) |
nonzero |
非零实数，(x\in\mathbb{R}-{0})。nonzero 数总是 real，且不能是 zero。 |
-> 实数``== (extended_nonzero & 有限) |
extended_nonzero |
扩展实数中不为零的成员，(x\in\overline{\mathbb{R}}-{0})。 | == (扩展实数 & !零) |
positive |
正实数，(x\in\mathbb{R}, x>0)。所有的 positive 数都是有限的，因此 oo 不是 positive。[正数] |
== (非负数 & 非零)``== (扩展正数 & 有限) |
nonnegative |
非负实数，(x\in\mathbb{R}, x\ge 0)。所有的 nonnegative 数都是有限的，因此 oo 不是 nonnegative。[正数] |
== (实数 & !负数) |
negative |
负实数，(x\in\mathbb{R}, x<0)。所有的 negative 数都是有限的，因此 -oo 不是 negative。[负数] |
== (非正数 & 非零)``== (扩展负数 & 有限) |
nonpositive |
非正实数，(x\in\mathbb{R}, x\le 0)。所有的 nonpositive 数都是有限的，因此 -oo 不是 nonpositive。[负数] |
== (实数 & !正数) |
extended_positive |
正的扩展实数，(x\in\overline{\mathbb{R}}, x>0)。extended_positive 数要么是 positive 要么是 oo。[扩展实数] |
== (扩展非负数 & 扩展非零) |
extended_nonnegative |
非负扩展实数，(x\in\overline{\mathbb{R}}, x\ge 0)。extended_nonnegative 数既可以是 nonnegative 也可以是 oo。[扩展实数] |
== (扩展实数 & !extended_negative) |
extended_negative |
负的扩展实数，(x\in\overline{\mathbb{R}}, x<0)。extended_negative 数要么是 negative 要么是 -oo。[扩展实数] |
== (扩展非正数 & 扩展非零) |
扩展非正数 |
非正扩展实数，(x\in\overline{\mathbb{R}}, x\le 0)。一个扩展非正数是非正数或-oo。[扩展实数] |
== (扩展实数 & !扩展正数) |
上述定义的参考资料
[交换性]
[无限] (1,2)
[反厄米] (1,2)
[复数]
[代数数]
[超越数]
[扩展实数] (1,2,3,4,5)
[实数]
[虚数]
[有理数]
[无理数]
[整数]
[奇偶性] (1,2)
zh.wikipedia.org/wiki/奇偶性 _(数学)
[质数]
[合数]
[零]
[正数] (1,2)
[负数] (1,2)
影响
假设系统使用推理规则来推断在创建符号时未明确指定的新谓词。
>>> x = Symbol('x', real=True, negative=False, zero=False)
>>> x.is_positive
True
虽然x没有明确声明为正数，但可以从明确给出的谓词推断出来。特别是推理规则之一是实数 == 负数 | 零 | 正数，所以如果实数为True且负数和零都为False，那么正数必须为True。
在实践中，假设推理规则意味着不必包含冗余的谓词，例如，正实数可以简单地声明为正：
>>> x1 = Symbol('x1', positive=True, real=True)
>>> x2 = Symbol('x2', positive=True)
>>> x1.is_real
True
>>> x2.is_real
True
>>> x1.assumptions0 == x2.assumptions0
True
合并不一致的谓词将导致错误：
>>> x = Symbol('x', commutative=False, real=True)
Traceback (most recent call last):
...
InconsistentAssumptions: {
algebraic: False,
commutative: False,
complex: False,
composite: False,
even: False,
extended_negative: False,
extended_nonnegative: False,
extended_nonpositive: False,
extended_nonzero: False,
extended_positive: False,
extended_real: False,
imaginary: False,
integer: False,
irrational: False,
negative: False,
noninteger: False,
nonnegative: False,
nonpositive: False,
nonzero: False,
odd: False,
positive: False,
prime: False,
rational: False,
real: False,
transcendental: False,
zero: False}, real=True
谓词的解释
尽管上表中定义了这些谓词，但值得花些时间考虑如何解释它们。首先，许多谓词名称所指的概念，如“零”，“素数”，“有理数”等在数学中有基本含义，但也可能有更广泛的含义。例如，在处理矩阵时，全零矩阵可能被称为“零”。假设系统中的谓词不允许这样的泛化。谓词 zero 严格保留给普通数 (0)。相反，矩阵具有一个 is_zero_matrix() 属性用于此目的（尽管该属性并不严格属于假设系统）：
>>> from sympy import Matrix
>>> M = Matrix([[0, 0], [0, 0]])
>>> M.is_zero
False
>>> M.is_zero_matrix
True
类似地，还有整数的泛化，如高斯整数，其对素数有不同的概念。假设系统中的 prime 谓词不包括这些，严格只指标准素数集合 (\mathbb{P} = {2, 3, 5, 7, 11, \cdots})。同样，integer 只意味着标准整数概念 (\mathbb{Z} = {0, \pm 1, \pm 2, \cdots})，rational 只意味着标准有理数概念 (\mathbb{Q})，依此类推。
这些谓词建立了类似于以复数开头的子集的体系，复数被认为是实数的超集，而实数又是有理数的超集，依此类推。子集的链条
[\mathbb{Z} \subset \mathbb{Q} \subset \mathbb{R} \subset \mathbb{C}]
对应于假设系统中的推理链条
integer -> rational -> real -> complex
一个“普通”的符号没有明确附加到这些集合中，甚至不知道它是有限的：
>>> x = Symbol('x')
>>> x.assumptions0
{'commutative': True}
>>> print(x.is_commutative)
True
>>> print(x.is_rational)
None
>>> print(x.is_complex)
None
>>> print(x.is_real)
None
>>> print(x.is_integer)
None
>>> print(x.is_finite)
None
对于 SymPy 来说，它很难知道如此一个甚至不知道是有限还是复数的符号可以做什么，因此通常最好明确地给符号一些假设。SymPy 的许多部分将会隐式地将这样的符号视为复数，在某些情况下，SymPy 将允许进行不严格的操作，尽管 x 不知道是有限的。从正式意义上讲，对于一个普通符号了解得很少，这使得涉及它的操作变得困难。
关于符号的定义可以产生很大的差异。例如，如果我们声明符号是整数，那么这意味着一系列其他谓词将有助于进一步的操作：
>>> n = Symbol('n', integer=True)
>>> n.assumptions0
{'algebraic': True,
'commutative': True,
'complex': True,
'extended_real': True,
'finite': True,
'hermitian': True,
'imaginary': False,
'infinite': False,
'integer': True,
'irrational': False,
'noninteger': False,
'rational': True,
'real': True,
'transcendental': False}
这些假设可以导致非常重要的简化，例如 integer=True 得到：
>>> from sympy import sin, pi
>>> n1 = Symbol('n1')
>>> n2 = Symbol('n2', integer=True)
>>> sin(n1 * pi)
sin(pi*n1)
>>> sin(n2 * pi)
0
用 (0) 替换整个表达式就像简化所能做到的一样好！
通常建议在任何符号上设置尽可能多的假设，以便尽可能简化表达式。一个常见的误解导致用 False 谓词定义一个符号，例如：
>>> x = Symbol('x', negative=False)
>>> print(x.is_negative)
False
>>> print(x.is_nonnegative)
None
>>> print(x.is_real)
None
>>> print(x.is_complex)
None
>>> print(x.is_finite)
None
如果意图是说 x 是一个不是正数的实数，那么需要明确说明。在已知符号是实数的情况下，谓词 positive=False 变得更有意义：
>>> x = Symbol('x', real=True, negative=False)
>>> print(x.is_negative)
False
>>> print(x.is_nonnegative)
True
>>> print(x.is_real)
True
>>> print(x.is_complex)
True
>>> print(x.is_finite)
True
声明为 Symbol('x', real=True, negative=False) 的符号等同于声明为 Symbol('x', nonnegative=True)。仅仅声明一个符号为 Symbol('x', positive=False) 并不能让假设系统对它有太多推断，因为普通符号并不被认为是有限的，甚至是复数。
关于 Symbol('x', complex=True) 和 Symbol('x', real=False) 存在一个相关的混淆。通常情况下，当其中任何一个被使用时，实际上不是我们想要的。首先要理解的是，所有的实数都是复数，因此用 real=True 创建的符号也会有 complex=True，而用 complex=True 创建的符号不会有 real=False。如果意图是创建一个既是复数又不是实数的复数，则应该使用 Symbol('x', complex=True, real=False)。另一方面，仅声明 real=False 是不足以推断 complex=True 的，因为知道它不是实数并不能告诉我们它是否有限或者是否是复数之外的完全不同的对象。
一个普通符号的定义是不知道它是否 finite 等，但是没有清晰的定义告诉我们它实际上应该代表什么。有时会想把它看作是一个“任意复数或可能是无穷大之一”，但是没有办法查询一个任意的（非符号）表达式以确定它是否符合这些条件。需要牢记的是，在 SymPy 代码库中以及可能在下游库中，还可以找到许多其他类型的数学对象，它们可能也有 commutative=True，但与普通数（在这种情况下，即使是 SymPy 的标准无穷大也被视为“普通”）完全不同。
唯一默认应用于符号的谓词是 commutative。我们还可以声明一个符号为 noncommutative，例如：
>>> x, y = symbols('x, y', commutative=False)
>>> z = Symbol('z') # defaults to commutative=True
>>> x*y + y*x
x*y + y*x
>>> x*z + z*x
2*z*x
注意这里，由于 x 和 y 都是非交换的，所以 x 和 y 不交换，即 x*y != y*x。另一方面，由于 z 是可交换的，所以 x 和 z 交换，即 x*z == z*x，尽管 x 是非交换的。
对于一个普通符号表示的解释不清楚，但是对于一个 commutative=False 的表达式的解释完全模糊。这样的表达式必然不是复数，也不是扩展实数或任何标准的无穷大（即使 zoo 是可交换的）。我们对这样一个表达式 代表 什么知之甚少。
其他 is_* 属性
SymPy 中有许多属性和特性，其名称以 is_ 开头，看起来类似于（旧）假设系统中使用的属性，但实际上并不属于假设系统。其中一些与假设系统的属性具有类似的含义和用法，例如上面显示的 is_zero_matrix() 属性。另一个例子是集合的 is_empty 属性：
>>> from sympy import FiniteSet, Intersection
>>> S1 = FiniteSet(1, 2)
>>> S1
{1, 2}
>>> print(S1.is_empty)
False
>>> S2 = Intersection(FiniteSet(1), FiniteSet(Symbol('x')))
>>> S2
Intersection({1}, {x})
>>> print(S2.is_empty)
None
is_empty 属性给出了一个模糊布尔值，指示一个 Set 是否为空集。在 S2 的示例中，不知道集合是否为空集，因为不知道 x 是否等于 1，所以 S2.is_empty 返回 None。对于集合，is_empty 属性起到类似于假设系统中 is_zero 属性的作用：is_empty 通常仅对 EmptySet 对象为 True，但仍然有助于区分 is_empty=False 和 is_empty=None 的情况。
尽管 is_zero_matrix 和 is_empty 用于与假设属性类似的目的，如 is_zero，但它们并不是（旧）假设系统的一部分。例如，没有关联的推理规则将 Set.is_empty 和 Set.is_finite_set 连接起来，因为推理规则是（旧）假设系统的一部分，该系统仅处理表格中列出的谓词。不可能声明一个具有例如 zero_matrix=False 的 MatrixSymbol，也没有 SetSymbol 类，但如果有的话，它不会有一个像 empty=False 这样的谓词理解系统。
属性 is_zero_matrix() 和 is_empty 类似于假设系统的属性，因为它们涉及表达式的语义方面。 还有许多其他属性，专注于结构方面，例如 is_Number，is_number()，is_comparable()。由于这些属性涉及表达式的结构方面，它们将始终返回 True 或 False，而不是具有可能为 None 的模糊布尔值。大写的属性如 is_Number 通常是 isinstance 检查的简写，例如：
>>> from sympy import Number, Rational
>>> x = Rational(1, 2)
>>> isinstance(x, Number)
True
>>> x.is_Number
True
>>> y = Symbol('y', rational=True)
>>> isinstance(y, Number)
False
>>> y.is_Number
False
Number 类是 Integer、Rational 和 Float 的超类，因此 Number 的任何实例代表具有已知值的具体数。像 y 这样用 rational=True 声明的符号可能代表与 x 相同的值，但它不是具有已知值的具体数，因此这是一种结构上而不是语义上的区别。例如 is_Number 属性有时在 SymPy 中用来取代 isinstance(obj, Number)，因为它们不会导致循环导入问题，并且检查 x.is_Number 可能比调用 isinstance 更快。
小写的 is_number 属性与 is_Number 非常不同。 is_number 属性对于任何可以通过 evalf() 数值评估为浮点复数的表达式都为 True：
>>> from sympy import I
>>> expr1 = I + sqrt(2)
>>> expr1
sqrt(2) + I
>>> expr1.is_number
True
>>> expr1.evalf()
1.4142135623731 + 1.0*I
>>> x = Symbol('x')
>>> expr2 = 1 + x
>>> expr2
x + 1
>>> expr2.is_number
False
>>> expr2.evalf()
x + 1.0
检查 expr.is_number 的主要原因是为了预测是否调用 evalf() 将完全评估。 is_comparable() 属性类似于 is_number()，但如果 is_comparable 返回 True，则表达式保证数值评估为 实数 Float。当 a.is_comparable 和 b.is_comparable 时，不等式 a < b 应该被解析为类似于 a.evalf() < b.evalf() 的形式。
在 SymPy 中，完整的 is_* 属性、属性和方法集合非常庞大。不过重要的是要清楚，只有那些在上面谓词表中列出的才是实际上的假设系统的一部分。只有那些涉及实现假设系统的 机制 的属性才会在下面解释。
实现假设处理程序
我们现在将通过实现一个 SymPy 符号函数的示例来说明如何查看内部使用的旧假设。SymPy 已经有一个对所有复数定义的 exp 函数，但我们将定义一个仅限于实数参数的 expreal 函数。
>>> from sympy import Function
>>> from sympy.core.logic import fuzzy_and, fuzzy_or
>>>
>>> class expreal(Function):
... """exponential function E**x restricted to the extended reals"""
...
... is_extended_nonnegative = True
...
... @classmethod
... def eval(cls, x):
... # Validate the argument
... if x.is_extended_real is False:
... raise ValueError("non-real argument to expreal")
... # Evaluate for special values
... if x.is_zero:
... return S.One
... elif x.is_infinite:
... if x.is_extended_negative:
... return S.Zero
... elif x.is_extended_positive:
... return S.Infinity
...
... @property
... def x(self):
... return self.args[0]
...
... def _eval_is_finite(self):
... return fuzzy_or([self.x.is_real, self.x.is_extended_nonpositive])
...
... def _eval_is_algebraic(self):
... if fuzzy_and([self.x.is_rational, self.x.is_nonzero]):
... return False
...
... def _eval_is_integer(self):
... if self.x.is_zero:
... return True
...
... def _eval_is_zero(self):
... return fuzzy_and([self.x.is_infinite, self.x.is_extended_negative])
Function.eval方法用于捕获函数的特殊值，以便我们可以在简化时返回不同的对象。当调用expreal(x)时，expreal.__new__类方法（在超类Function中定义）将调用expreal.eval(x)。如果expreal.eval返回的不是None，则将返回该值，而不是未求值的expreal(x)：
>>> from sympy import oo
>>> expreal(1)
expreal(1)
>>> expreal(0)
1
>>> expreal(-oo)
0
>>> expreal(oo)
oo
注意，expreal.eval方法不使用==比较参数。特殊值是使用假设系统来查询参数的属性。这意味着expreal方法也可以对具有匹配属性的不同形式的表达式进行评估，例如：
>>> x = Symbol('x', extended_negative=True, infinite=True)
>>> x
x
>>> expreal(x)
0
当然，假设系统只能解析有限数量的特殊值，因此大多数eval方法也会使用==检查一些特殊值，但最好检查例如x.is_zero而不是x==0。
还要注意，expreal.eval方法验证参数是否为实数。我们希望允许(\pm\infty)作为expreal的参数，因此我们检查extended_real而不是real。如果参数不是扩展实数，则会引发错误：
>>> expreal(I)
Traceback (most recent call last):
...
ValueError: non-real argument to expreal
重要的是，我们检查x.is_extended_real is False而不是not x.is_extended_real，这意味着我们只在参数绝对不是扩展实数时拒绝它：如果x.is_extended_real返回None，则不会拒绝参数。允许x.is_extended_real=None的第一个原因是可以将普通符号用于expreal。第二个原因是，在即使参数明确为实数的情况下，假设查询也可能始终返回None，例如：
>>> x = Symbol('x')
>>> print(x.is_extended_real)
None
>>> expreal(x)
expreal(x)
>>> expr = (1 + I)/sqrt(2) + (1 - I)/sqrt(2)
>>> print(expr.is_extended_real)
None
>>> expr.expand()
sqrt(2)
>>> expr.expand().is_extended_real
True
>>> expreal(expr)
expreal(sqrt(2)*(1 - I)/2 + sqrt(2)*(1 + I)/2)
在expreal.eval中验证参数确实意味着在传递evaluate=False时不会验证参数，但实际上没有更好的位置来执行验证：
>>> expreal(I, evaluate=False)
expreal(I)
expreal类的extended_nonnegative类属性和expreal类上的_eval_is_*方法实现了对expreal实例在假设系统中的查询：
>>> expreal(2)
expreal(2)
>>> expreal(2).is_finite
True
>>> expreal(2).is_integer
False
>>> expreal(2).is_rational
False
>>> expreal(2).is_algebraic
False
>>> z = expreal(-oo, evaluate=False)
>>> z
expreal(-oo)
>>> z.is_integer
True
>>> x = Symbol('x', real=True)
>>> expreal(x)
expreal(x)
>>> expreal(x).is_nonnegative
True
假设系统使用相应的处理程序expreal._eval_is_finite解析像expreal(2).is_finite这样的查询，还使用蕴含规则。例如，已知expreal(2).is_rational为False，因为expreal(2)._eval_is_algebraic返回False，并且存在一个蕴含规则rational -> algebraic。这意味着在这种情况下，可以通过_eval_is_algebraic处理程序解析is_rational查询。实际上，最好不要为每个可能的谓词实现假设处理程序，而是尝试识别一组最小的处理程序，可以尽可能少地检查以解析尽可能多的查询：
另一个需要注意的地方是，_eval_is_* 方法仅对参数 x 进行假设查询，并不对 self 进行任何假设查询。在同一对象上递归进行假设查询会干扰假设蕴含解析器，可能导致非确定性行为，因此不应使用它们（SymPy 代码库中存在这样的例子，但应予以移除）。
许多 expreal 方法隐式返回 None。这是假设系统中的常见模式。eval 方法和 _eval_is_* 方法都可以返回 None，并且通常会返回 None。Python 函数如果在没有达到 return 语句的情况下结束，会隐式返回 None。我们通过省略许多 if 语句的 else 子句，并允许隐式返回 None 来利用这一点。在跟踪这些方法的控制流时，重要的是要记住，任何查询的属性都可能返回 True、False 或 None，并且如果所有条件失败，任何函数都会隐式返回 None。
假设系统的机制
注意
本节描述了在未来 SymPy 版本中可能会发生变化的内部细节。
本节将解释假设系统的内部工作原理。重要的是要理解，这些内部工作原理是实现细节，可能会从一个 SymPy 版本变化到另一个版本。此解释适用于 SymPy 1.7 版本。尽管（旧）假设系统存在许多限制（在下一节讨论），但它是一个成熟的系统，在 SymPy 中被广泛使用，并且已经针对当前用途进行了优化。假设系统在大多数 SymPy 操作中隐式使用，以控制基本表达式的评估。
在 SymPy 进程中，假设系统的实现经历了几个阶段，最终导致假设系统中单个查询的评估。简而言之，这些阶段包括：
-
在导入时，
sympy/core/assumptions.py中定义的假设规则被处理成一个规范形式，以便高效地应用蕴含规则。这在 SymPy 被导入时仅发生一次，甚至在Basic类被定义之前。 -
Basic.__init_subclass__方法将后处理每个Basic子类，以添加所需的属性，用于假设查询。这也向类添加了default_assumptions属性。每次定义Basic子类时（导入其所在的模块时），都会发生这种情况。 -
每个
Basic实例最初使用default_assumptions类属性。当对Basic实例进行假设查询时，在第一次实例中，查询将从类的default_assumptions中得到答案。 -
如果类的
default_assumptions中对于假设查询没有缓存值，则将复制默认假设以创建实例的假设缓存。然后调用_ask()函数来解析查询，该函数首先将调用相关实例处理程序_eval_is方法。如果处理程序返回非None，则结果将被缓存并返回。 -
如果处理程序不存在或返回
None，则将尝试使用推理解析器。这将以随机顺序枚举所有可能用于根据推理规则解析查询的谓词组合。在每种情况下，将调用处理程序的_eval_is方法以查看是否返回非None。如果任何处理程序和推理规则的组合导致查询的确定结果，则将该结果缓存到实例缓存中并返回。 -
最后，如果推理解析器未能解析查询，则认为查询无法解决。查询的
None值将被缓存在实例缓存中并返回。
sympy/core/assumptions.py 中定义的假设规则以 real == negative | zero | positive 的形式给出。当导入这个模块时，这些规则会被转换成一个名为 _assume_rules 的 FactRules 实例。这将预处理含义规则，转换成可以用于推理解析的 “A” 和 “B” 规则形式。这在 sympy/core/facts.py 中的代码中有详细说明。我们可以像这样直接访问这个内部对象（完整输出略）：
>>> from sympy.core.assumptions import _assume_rules
>>> _assume_rules.defined_facts
{'algebraic',
'antihermitian',
'commutative',
'complex',
'composite',
'even',
...
>>> _assume_rules.full_implications
defaultdict(set,
{('extended_positive', False): {('composite', False),
('positive', False),
('prime', False)},
('finite', False): {('algebraic', False),
('complex', False),
('composite', False),
...
Basic.__init_subclass__ 方法将检查每个 Basic 类的属性，看看是否定义了任何与假设相关的属性。这些属性的一个示例是在 expreal 类中定义的 is_extended_nonnegative = True 属性。这些属性的含义将用于预先计算任何静态可知的假设。例如，is_extended_nonnegative=True 暗示 real=True 等。为类创建一个 StdFactKB 实例，该实例存储这些在此阶段已知的假设值。将 StdFactKB 实例分配为类属性 default_assumptions。我们可以通过以下方式看到这一点：
>>> from sympy import Expr
...
>>> class A(Expr):
... is_positive = True
...
... def _eval_is_rational(self):
... # Let's print something to see when this method is called...
... print('!!! calling _eval_is_rational')
... return True
...
>>> A.is_positive
True
>>> A.is_real # inferred from is_positive
True
尽管在类 A 中仅定义了 is_positive，它还具有诸如 is_real 等从 is_positive 推断而来的属性。类 A 的所有这类假设集合可以在 default_assumptions 中看到，它看起来像一个 dict，但实际上是一个 StdFactKB 实例：
>>> type(A.default_assumptions)
<class 'sympy.core.assumptions.StdFactKB'>
>>> A.default_assumptions
{'commutative': True,
'complex': True,
'extended_negative': False,
'extended_nonnegative': True,
'extended_nonpositive': False,
'extended_nonzero': True,
'extended_positive': True,
'extended_real': True,
'finite': True,
'hermitian': True,
'imaginary': False,
'infinite': False,
'negative': False,
'nonnegative': True,
'nonpositive': False,
'nonzero': True,
'positive': True,
'real': True,
'zero': False}
当创建任何Basic子类的实例时，Basic.__new__将分配其 _assumptions 属性，该属性最初将是对 cls.default_assumptions 的引用，这在同一类的所有实例之间是共享的。实例将使用此属性来解析任何假设查询，直到无法给出明确结果为止，此时将创建 cls.default_assumptions 的副本，并分配给实例的 _assumptions 属性。该副本将用作缓存，用于存储由其 _eval_is 处理程序计算的实例的任何结果。
当 _assumptions 属性未能提供相关结果时，是调用 _eval_is 处理程序的时候了。此时会调用 _ask() 函数。_ask() 函数将首先尝试通过调用相应的方法，如 _eval_is_rational 来解析诸如 is_rational 的查询。如果返回非空，则结果将存储在 _assumptions 中，并计算并存储该结果的任何推论。此时，查询将被解析并返回值。
>>> a = A()
>>> a._assumptions is A.default_assumptions
True
>>> a.is_rational
!!! calling _eval_is_rational
True
>>> a._assumptions is A.default_assumptions
False
>>> a._assumptions # rational now shows as True
{'algebraic': True,
'commutative': True,
'complex': True,
'extended_negative': False,
'extended_nonnegative': True,
'extended_nonpositive': False,
'extended_nonzero': True,
'extended_positive': True,
'extended_real': True,
'finite': True,
'hermitian': True,
'imaginary': False,
'infinite': False,
'irrational': False,
'negative': False,
'nonnegative': True,
'nonpositive': False,
'nonzero': True,
'positive': True,
'rational': True,
'real': True,
'transcendental': False,
'zero': False}
如果例如 _eval_is_rational 不存在或返回 None，那么 _ask() 将尝试所有可能性来使用推论规则和任何其他处理程序方法，如 _eval_is_integer、_eval_is_algebraic 等，这些方法可能能够对原始查询给出答案。如果任何方法导致已知原始查询的确定结果，则返回该结果。否则，一旦用于解析查询的处理程序和推论规则的所有可能性都耗尽，将缓存并返回 None。
>>> b = A()
>>> b.is_algebraic # called _eval_is_rational indirectly
!!! calling _eval_is_rational
True
>>> c = A()
>>> print(c.is_prime) # called _eval_is_rational indirectly
!!! calling _eval_is_rational
None
>>> c._assumptions # prime now shows as None
{'algebraic': True,
'commutative': True,
'complex': True,
'extended_negative': False,
'extended_nonnegative': True,
'extended_nonpositive': False,
'extended_nonzero': True,
'extended_positive': True,
'extended_real': True,
'finite': True,
'hermitian': True,
'imaginary': False,
'infinite': False,
'irrational': False,
'negative': False,
'nonnegative': True,
'nonpositive': False,
'nonzero': True,
'positive': True,
'prime': None,
'rational': True,
'real': True,
'transcendental': False,
'zero': False}
注意
在 _ask() 函数中，处理程序的调用顺序是随机的，这意味着此时的执行是非确定性的。只要所有不同的处理程序方法保持一致（即没有错误），最终结果仍将是确定性的。然而，如果存在两个处理程序不一致的 bug，则可能导致非确定性行为，因为此随机化可能导致在多次运行同一程序时以不同的顺序调用处理程序。
限制
合并带有 or 的谓词
在旧的假设中，我们可以轻松地在创建符号时将谓词与and组合，例如：
>>> x = Symbol('x', integer=True, positive=True)
>>> x.is_positive
True
>>> x.is_integer
True
我们还可以轻松地查询两个条件是否共同满足
>>> fuzzy_and([x.is_positive, x.is_integer])
True
>>> x.is_positive and x.is_integer
True
然而，在旧的假设中，无法创建具有或组合的Symbol。例如，如果我们想要表达“x 是正数或 x 是整数”，则无法创建具有这些假设的Symbol。
也不可能根据或进行假设查询，例如“expr 是否为正数或整数表达式”。我们可以使用例如
>>> fuzzy_or([x.is_positive, x.is_integer])
True
但是，如果对x的所有了解只是它可能是正数或者是一个负整数，那么x.is_positive和x.is_integer这两个查询都会返回None。这意味着查询变成了
>>> fuzzy_or([None, None])
然后也会返回None。
不同符号之间的关系
旧假设系统的一个基本限制是，所有显式假设都是一个单独符号的属性。在这个系统中，没有办法假设两个符号之间的关系。最常见的请求之一是能够假设类似x < y这样的内容，但在旧的假设中甚至没有办法指定这一点。
新的假设具有理论上的能力，可以指定关系性假设。然而，利用该信息的算法尚未实现，并且还未决定指定关系性假设的确切 API。
符号与模糊布尔
本页描述了 SymPy 中符号化 Boolean 的含义，以及它与用于 SymPy 许多部分的三值模糊布尔的关系。还讨论了在编写使用三值逻辑的代码时可能出现的一些常见问题，以及如何正确处理它们。
符号布尔与三值布尔
像 x.ispositive 这样的假设查询会给出三值模糊布尔的 True、False 或 None 结果 [1]。这些是低级别的 Python 对象，而不是 SymPy 的符号化 Boolean 表达式。
>>> from sympy import Symbol, symbols
>>> xpos = Symbol('xpos', positive=True)
>>> xneg = Symbol('xneg', negative=True)
>>> x = Symbol('x')
>>> print(xpos.is_positive)
True
>>> print(xneg.is_positive)
False
>>> print(x.is_positive)
None
作为三值模糊布尔的 None 结果应该解释为“可能”或“未知”。
在 SymPy 中，可以在不等式中找到一个符号 Boolean 类。当不等式不能确定为真或假时，Boolean 可以以符号方式表示不确定的结果：
>>> xpos > 0
True
>>> xneg > 0
False
>>> x > 0
x > 0
>>> type(x > 0)
<class 'sympy.core.relational.StrictGreaterThan'>
最后一个例子显示了当不等式不确定时会发生什么：我们得到一个 StrictGreaterThan 的实例，它表示不等式作为一个符号表达式。内部当尝试评估像 a > b 这样的不等式时，SymPy 将计算 (a - b).is_extended_positive。如果结果是 True 或 False，则会返回 SymPy 的符号 S.true 或 S.false。如果结果是 None，则会返回一个未求值的 StrictGreaterThan，如上面的 x > 0 所示。
并不明显，像 xpos > 0 这样的查询返回的是 S.true 而不是 True，因为这两个对象在显示上相同，但我们可以使用 Python 的 is 运算符来验证：
>>> from sympy import S
>>> xpos.is_positive is True
True
>>> xpos.is_positive is S.true
False
>>> (xpos > 0) is True
False
>>> (xpos > 0) is S.true
True
在 SymPy 中，没有 None 的通用符号类似物。在低级别的假设查询中返回 None 的情况下，符号查询将导致一个未求值的符号化 Boolean （例如，x > 0）。我们可以将符号 Boolean 用作符号表达式的一部分，如 Piecewise：
>>> from sympy import Piecewise
>>> p = Piecewise((1, x > 0), (2, True))
>>> p
Piecewise((1, x > 0), (2, True))
>>> p.subs(x, 3)
1
这里 p 表示一个表达式，如果 x > 0 则等于 1，否则将等于 2。未评估的 Boolean 不等式 x > 0 表示决定符号表达式值的条件。当我们为 x 替换一个值时，不等式将解析为 S.true，然后 Piecewise 可以评估为 1 或 2。
当使用模糊布尔值而不是符号 Boolean 时，同样不起作用：
>>> p2 = Piecewise((1, x.is_positive), (2, True))
Traceback (most recent call last):
...
TypeError: Second argument must be a Boolean, not `NoneType`
Piecewise 不能将 None 用作条件，因为与不等式 x > 0 不同，它没有提供任何信息。使用不等式时可以在将来决定条件是否可能为 True 或 False，一旦知道 x 的值。None 的值不能以这种方式使用，因此被拒绝。
注意
我们可以在 Piecewise 中使用 True，因为 True sympifies 为 S.true。将 None sympify 只会再次得到 None，这不是一个有效的符号 SymPy 对象。
SymPy 中有许多其他符号 Boolean 类型。关于模糊布尔值和符号 Boolean 之间的差异的同样考虑适用于所有其他 SymPy Boolean 类型。举一个不同的例子，有 Contains，它表示对象包含在集合中的陈述：
>>> from sympy import Reals, Contains
>>> x = Symbol('x', real=True)
>>> y = Symbol('y')
>>> Contains(x, Reals)
True
>>> Contains(y, Reals)
Contains(y, Reals)
>>> Contains(y, Reals).subs(y, 1)
True
对应于 Contains 的 Python 操作符是 in。in 的一个怪异之处在于它只能评估为 bool（True 或 False），所以如果结果是不确定的，则会引发异常：
>>> from sympy import I
>>> 2 in Reals
True
>>> I in Reals
False
>>> x in Reals
True
>>> y in Reals
Traceback (most recent call last):
...
TypeError: did not evaluate to a bool: (-oo < y) & (y < oo)
可以通过使用 Contains(x, Reals) 或 Reals.contains(x) 而不是 x in Reals 来避免异常。
使用模糊布尔值的三值逻辑
无论我们使用模糊布尔值还是符号 Boolean，我们始终需要意识到查询可能是不确定的。如何编写处理此问题的代码在两种情况下是不同的。我们先看看模糊布尔值。
考虑以下函数：
>>> def both_positive(a, b):
... """ask whether a and b are both positive"""
... if a.is_positive and b.is_positive:
... return True
... else:
... return False
both_positive函数应告诉我们a和b是否都是正数。然而，如果任何一个is_positive查询返回None，both_positive函数将失败：
>>> print(both_positive(S(1), S(1)))
True
>>> print(both_positive(S(1), S(-1)))
False
>>> print(both_positive(S(-1), S(-1)))
False
>>> x = Symbol('x') # may or may not be positive
>>> print(both_positive(S(1), x))
False
注意
我们需要使用S对这个函数的参数进行符号化，因为假设仅在 SymPy 对象上定义，而不是普通的 Python int对象上定义。
在这里，False是不正确的，因为可能x是正数，此时两个参数都可能是正数。我们得到False是因为x.is_positive返回None，Python 会将None视为“假值”。
为了正确处理所有可能情况，我们需要分开识别True和False情况的逻辑。一个改进的函数可能是：
>>> def both_positive_better(a, b):
... """ask whether a and b are both positive"""
... if a.is_positive is False or b.is_positive is False:
... return False
... elif a.is_positive is True and b.is_positive is True:
... return True
... else:
... return None
这个函数现在可以处理所有情况的True、False或None，对于a和b都是如此，并且总是返回一个模糊布尔值，表示语句“a和b都是正数”是真、假还是未知：
>>> print(both_positive_better(S(1), S(1)))
True
>>> print(both_positive_better(S(1), S(-1)))
False
>>> x = Symbol('x')
>>> y = Symbol('y', positive=True)
>>> print(both_positive_better(S(1), x))
None
>>> print(both_positive_better(S(-1), x))
False
>>> print(both_positive_better(S(1), y))
True
使用模糊布尔值时需要小心的另一种情况是与 Python 的not运算符的否定，例如：
>>> x = Symbol('x')
>>> print(x.is_positive)
None
>>> not x.is_positive
True
模糊布尔值None的正确否定再次是None。如果我们不知道语句“x是正数”是True还是False，那么我们也不知道其否定“x不是正数”是True还是False。之所以会得到True，是因为None被视为“假值”。当None与诸如not这样的逻辑运算符一起使用时，它首先会转换为bool，然后再取反：
>>> bool(None)
False
>>> not bool(None)
True
>>> not None
True
当None被视为假值时，如果使用正确，这一点是有用的。例如，如果我们只想在x被认为是正数的情况下执行某些操作，我们可以这样做：
>>> x = Symbol('x', positive=True)
>>> if x.is_positive:
... print("x is definitely positive")
... else:
... print("x may or may not be positive")
x is definitely positive
只要我们理解备用条件分支指的是两种情况（False和None），这种写法可以是编写条件语句的一个有用方式。当我们确实需要区分所有情况时，我们需要使用诸如x.is_positive is False之类的语句。但需要注意的是，当使用 Python 的二进制逻辑运算符如not或and处理模糊布尔值时，它们不能正确处理不定情况。
实际上，SymPy 有内部函数专门设计用于正确处理模糊布尔值：
>>> from sympy.core.logic import fuzzy_not, fuzzy_and
>>> print(fuzzy_not(True))
False
>>> print(fuzzy_not(False))
True
>>> print(fuzzy_not(None))
None
>>> print(fuzzy_and([True, True]))
True
>>> print(fuzzy_and([True, None]))
None
>>> print(fuzzy_and([False, None]))
False
使用fuzzy_and函数，我们可以更简单地编写both_positive函数：
>>> def both_positive_best(a, b):
... """ask whether a and b are both positive"""
... return fuzzy_and([a.is_positive, b.is_positive])
利用fuzzy_and、fuzzy_or和fuzzy_not编写更简洁的代码不仅可以减少逻辑错误的可能性，还可以使代码看起来更像普通二进制逻辑的情况。
三值逻辑与符号布尔值
在使用符号Boolean而不是模糊布尔值时，None被隐式视为假值的问题不会出现，因此如果不小心处理，不会出现逻辑错误。然而，代之以不定情况通常会导致异常被抛出，如果不小心处理的话。
这次我们将尝试使用符号 Boolean 实现 both_positive 函数：
>>> def both_positive(a, b):
... """ask whether a and b are both positive"""
... if a > 0 and b > 0:
... return S.true
... else:
... return S.false
第一个区别是，我们返回符号 Boolean 对象 S.true 和 S.false 而不是 True 和 False。第二个区别是，我们测试例如 a > 0 而不是 a.is_positive。尝试这样做我们得到
>>> both_positive(1, 2)
True
>>> both_positive(-1, 1)
False
>>> x = Symbol('x') # may or may not be positive
>>> both_positive(x, 1)
Traceback (most recent call last):
...
TypeError: cannot determine truth value of Relational
现在的情况是，当 x 不知道是否为正数或非正数时，测试 x > 0 会引发异常。更准确地说，x > 0 不会引发异常，但是 if x > 0 会，这是因为 if 语句隐式调用 bool(x > 0)，而后者会引发异常。
>>> x > 0
x > 0
>>> bool(x > 0)
Traceback (most recent call last):
...
TypeError: cannot determine truth value of Relational
>>> if x > 0:
... print("x is positive")
Traceback (most recent call last):
...
TypeError: cannot determine truth value of Relational
Python 表达式 x > 0 创建了一个 SymPy Boolean。因为在这种情况下，Boolean 不能评估为 True 或 False，所以我们得到了一个未评估的 StrictGreaterThan。试图强制转换为 bool 类型，即 bool(x > 0) 会引发异常。这是因为普通的 Python bool 类型必须是 True 或 False，而在这种情况下，这两者都不确定是正确的。
当使用符号 Boolean 与 and、or 或 not 时会出现相同类型的问题。解决方案是使用 SymPy 的符号 And、Or 和 Not 或者等效的 Python 位逻辑运算符 &、| 和 ~：
>>> from sympy import And, Or, Not
>>> x > 0
x > 0
>>> x > 0 and x < 1
Traceback (most recent call last):
...
TypeError: cannot determine truth value of Relational
>>> And(x > 0, x < 1)
(x > 0) & (x < 1)
>>> (x > 0) & (x < 1)
(x > 0) & (x < 1)
>>> Or(x < 0, x > 1)
(x > 1) | (x < 0)
>>> Not(x < 0)
x >= 0
>>> ~(x < 0)
x >= 0
如前所述，如果避免直接在 if、and、or 或 not 中使用 SymPy 的符号 Boolean，我们可以创建一个更好的 both_positive 版本，而是测试 Boolean 是否评估为 S.true 或 S.false：
>>> def both_positive_better(a, b):
... """ask whether a and b are both positive"""
... if (a > 0) is S.false or (b > 0) is S.false:
... return S.false
... elif (a > 0) is S.true and (b > 0) is S.true:
... return S.true
... else:
... return And(a > 0, b > 0)
现在，使用这个版本，我们不会得到任何异常，如果结果是不确定的，我们将得到一个表示语句“a 和 b 都是正数”的符号 Boolean：
>>> both_positive_better(S(1), S(2))
True
>>> both_positive_better(S(1), S(-1))
False
>>> x, y = symbols("x, y")
>>> both_positive_better(x, y + 1)
(x > 0) & (y + 1 > 0)
>>> both_positive_better(x, S(3))
x > 0
最后一个案例显示，实际上使用已知为真的条件与 And 简化了 And。实际上我们有
>>> And(x > 0, 3 > 0)
x > 0
>>> And(4 > 0, 3 > 0)
True
>>> And(-1 > 0, 3 > 0)
False
这意味着我们可以改进both_positive_better。所有不同情况都是不必要的。我们可以简单地返回And，并让它在可能的情况下简化：
>>> def both_positive_best(a, b):
... """ask whether a and b are both positive"""
... return And(a > 0, b > 0)
现在这将适用于任何符号实对象，并产生一个符号结果。我们还可以替换结果，看看它如何适用于特定值：
>>> both_positive_best(2, 1)
True
>>> both_positive_best(-1, 2)
False
>>> both_positive_best(x, 3)
x > 0
>>> condition = both_positive_best(x/y, x + y)
>>> condition
(x + y > 0) & (x/y > 0)
>>> condition.subs(x, 1)
(1/y > 0) & (y + 1 > 0)
>>> condition.subs(x, 1).subs(y, 2)
True
当处理符号Boolean对象时，最好尽可能避免尝试使用if/else和其他逻辑运算符如and等进行分支。相反，考虑计算条件并将其作为变量传递。基本的符号操作如And、Or和Not可以为您处理逻辑。
脚注
编写自定义函数
本指南将描述如何在 SymPy 中创建自定义函数类。自定义用户定义函数使用与 SymPy 中包含的函数相同的机制，如函数中包含的常见初等函数，例如exp()或sin()，特殊函数如gamma()或Si()，以及组合函数和数论函数，如factorial()或primepi()。因此，本指南既是为希望定义自己自定义函数的最终用户提供指南，也是为希望扩展 SymPy 中包含的函数的 SymPy 开发人员提供指南。
本指南描述了如何定义复值函数，即将(\mathbb{C}^n)的子集映射到(\mathbb{C})的函数。接受或返回复数以外对象的函数应该是另一个类的子类，比如Boolean、MatrixExpr、Expr或Basic。这里写的一些内容适用于一般的Basic或Expr子类，但其中大部分仅适用于Function子类。
简单情况：完全符号化或完全评估
在深入研究自定义函数的更高级功能之前，我们应该提到两种常见情况，一个是函数完全符号化的情况，另一个是函数完全评估的情况。这两种情况都有比本指南中描述的完整机制更简单的替代方法。
完全符号情况
如果您的函数f没有您想要在其上定义的数学属性，并且不应在任何参数上进行评估，则可以使用Function('f')创建一个未定义的函数
>>> from sympy import symbols, Function
>>> x = symbols('x')
>>> f = Function('f')
>>> f(x)
f(x)
>>> f(0)
f(0)
这在解决 ODEs 时非常有用。
如果您只希望创建一个仅用于不同化目的依赖于另一个符号的符号，则这也是有用的。默认情况下，SymPy 假设所有符号彼此独立：
>>> from sympy.abc import x, y
>>> y.diff(x)
0
要创建一个依赖于另一个符号的符号，您可以使用明确依赖于该符号的函数。
>>> y = Function('y')
>>> y(x).diff(x)
Derivative(y(x), x)
如果您希望函数具有其他行为，例如具有自定义导数或在某些参数上进行评估，则应创建一个自定义Function子类，如下所述。但是，未定义的函数确实支持一个附加功能，即可以使用与符号相同的语法来定义它们的假设。这定义了函数输出的假设，而不是输入（即定义了函数的范围，而不是定义其域）。
>>> g = Function('g', real=True)
>>> g(x)
g(x)
>>> g(x).is_real
True
要使函数的假设依赖于其输入方式，您应创建一个自定义的Function子类，并如下所述定义假设处理程序。 ### 完全评估情况
在另一端的函数是始终评估为某些内容的函数，无论其输入如何。这些函数从不以未评估的符号形式如f(x)留下。
在这种情况下，您应该使用使用def关键字创建一个普通的 Python 函数：
>>> def f(x):
... if x == 0:
... return 0
... else:
... return x + 1
>>> f(0)
0
>>> f(1)
2
>>> f(x)
x + 1
如果您发现自己在Function子类上定义了一个eval()方法，其中您总是返回一个值，而不是返回None，那么考虑只是使用普通的 Python 函数，因为在这种情况下使用符号Function子类没有任何好处（参见下面的 eval()最佳实践部分）
注意，在许多情况下，这些函数可以直接使用 SymPy 类表示。例如，上述函数可以使用Piecewise进行符号表示。可以使用subs()对Piecewise表达式进行特定x值的评估。
>>> from sympy import Piecewise, Eq, pprint
>>> f = Piecewise((0, Eq(x, 0)), (x + 1, True))
>>> pprint(f, use_unicode=True)
⎧ 0 for x = 0
⎨
⎩x + 1 otherwise
>>> f.subs(x, 0)
0
>>> f.subs(x, 1)
2
像 Piecewise 这样的完全符号表示具有准确表示符号值的优势。例如，在上述 Python 的 def 定义 f 中，f(x) 隐式地假定 x 是非零的。Piecewise 版本会正确处理这种情况，并且不会在 x 不为零时评估到 (x \neq 0) 的情况。
如果您希望函数不仅进行评估，而且总是评估为数值，还有另一种选择，那就是使用lambdify()。这将把 SymPy 表达式转换为可以使用 NumPy 进行评估的函数。
>>> from sympy import lambdify
>>> func = lambdify(x, Piecewise((0, Eq(x, 0)), (x + 1, True)))
>>> import numpy as np
>>> func(np.arange(5))
array([0., 2., 3., 4., 5.])
最终，选择正确的工具取决于您要做什么以及您想要的确切行为。 ## 创建自定义函数
创建自定义函数的第一步是子类化Function。子类的名称将是函数的名称。然后，根据您想要提供的功能，应该在这个子类上定义不同的方法。
作为本文档的一个激励性例子，让我们创建一个表示versine 函数的自定义函数类。Versine 是一个三角函数，历史上与更熟悉的正弦和余弦函数一起使用。今天很少使用。Versine 可以通过下面的恒等式来定义
[\operatorname{versin}(x) = 1 - \cos(x).]
SymPy 不包括 versine，因为它在现代数学中很少使用，而且可以很容易地用更熟悉的余弦来定义。
让我们从子类化 Function 开始。
>>> class versin(Function):
... pass
此时，versin 没有定义任何行为。它与我们上面讨论过的未定义函数非常相似。请注意，versin 是一个类，versin(x) 是这个类的一个实例。
>>> versin(x)
versin(x)
>>> isinstance(versin(x), versin)
True
注意
下面描述的所有方法都是可选的。如果您希望定义特定的行为，可以包含它们，但如果省略它们，SymPy 将默认保持未评估状态。例如，如果您不定义微分，diff() 将只返回一个未评估的Derivative。
使用 eval() 定义自动评估
我们可能希望在自定义函数上定义的第一件事情是自动评估，即在返回实际值而不是保持未评估状态时的情况。
这是通过定义类方法eval()完成的。eval()应该接受函数的参数并返回一个值或None。如果返回None，则函数在那种情况下将保持未评估状态。这也有助于定义函数的签名（默认情况下，没有eval()方法，Function子类将接受任意数量的参数）。
对于我们的函数versin，我们可能会回忆起当整数n时，(\cos(n\pi) = (-1)^n)，因此(\operatorname{versin}(n\pi) = 1 - (-1)^n.) 当传递整数倍的pi时，我们可以使versin自动评估为这个值：
>>> from sympy import pi, Integer
>>> class versin(Function):
... @classmethod
... def eval(cls, x):
... # If x is an integer multiple of pi, x/pi will cancel and be an Integer
... n = x/pi
... if isinstance(n, Integer):
... return 1 - (-1)**n
>>> versin(pi)
2
>>> versin(2*pi)
0
在这里，我们利用了 Python 函数如果没有显式返回值，则自动返回None的事实。因此，在未触发if isinstance(n, Integer)语句的情况下，eval()返回None，并且versin保持未评估状态。
>>> versin(x*pi)
versin(pi*x)
注意
Function子类不应重新定义__new__或__init__。如果要实现eval()无法实现的行为，可能更合理的是子类化Expr而不是Function。
eval()可以接受任意数量的参数，包括带有*args和可选关键字参数的任意数量。函数的.args始终是用户传入的参数。例如
>>> class f(Function):
... @classmethod
... def eval(cls, x, y=1, *args):
... return None
>>> f(1).args
(1,)
>>> f(1, 2).args
(1, 2)
>>> f(1, 2, 3).args
(1, 2, 3)
最后，请注意，一旦定义了evalf()，浮点输入的自动评估就会自动发生，因此你不需要在eval()中显式处理它。
eval()的最佳实践
在定义eval()方法时，存在一些常见的反模式，应该避免。
-
不要只返回表达式。
在上面的例子中，我们可能会被诱惑写
>>> from sympy import cos >>> class versin(Function): ... @classmethod ... def eval(cls, x): ... # !! Not actually a good eval() method !! ... return 1 - cos(x)然而，这将导致
versin(x)始终返回1 - cos(x)，无论x是什么。如果你只想要一个快捷方式到1 - cos(x)，那没问题，但是更简单和更明确的方法是像上面描述的使用 Python 函数。如果我们像这样定义versin，它实际上永远不会表示为versin(x)，并且我们在versin类下面定义的任何其他行为都不会起作用，因为只有当返回的对象实际上是versin实例时，我们定义的其他行为才适用。例如，versin(x).diff(x)实际上只是(1 - cos(x)).diff(x)，而不是调用我们在下面定义的 fdiff()方法。关键点
eval()的目的不是数学上定义函数是什么，而是指定在哪些输入下它应该自动评估。 函数的数学定义是通过下面概述的各种数学属性的规范来确定的，比如 numerical evaluation，differentiation 等方法。如果你发现自己在这样做，请考虑你实际想要达到的目标。如果你只想为一个表达式定义一个简短的函数，最简单的方法就是定义一个 Python 函数。如果你真的想要一个符号函数，想一想你希望它在什么时候评估为其他值，以及什么时候保持不变。一种选择是在
eval()中使你的函数保持未评估状态，并定义一个doit()方法来评估它。 -
避免过多的自动评估。
建议最小化
eval()自动评估的内容。通常最好将更高级的简化放在其他方法中，如doit()。记住，无论你为自动评估定义什么，它都总是会进行评估。[1] 如前一点所述，如果你评估每个值，那么首先拥有符号函数就没有多大意义。例如，我们可能会试图在eval()中对versin进行一些三角恒等式的评估，但这些恒等式将始终被评估，并且无法表示恒等式的一半。还应避免在
eval()中执行计算速度慢的操作。SymPy 通常假设创建表达式是廉价的，如果不是这样，可能会导致性能问题。最后，建议避免根据假设在
eval()中进行自动评估。相反，eval()通常只评估显式的数值特定值，并对其他情况返回None。你可能已经注意到在上面的例子中我们使用了isinstance(n, Integer)而不是使用假设系统检查n.is_integer。我们本可以这样做，这样versin(n*pi)会被评估，即使n = Symbol('n', integer=True)。但这是一个情况，我们可能并不总是希望发生评估，如果n是一个更复杂的表达式，使用n.is_integer可能计算代价更高。考虑一个例子。使用恒等式 (\cos(x + y) = \cos(x)\cos(y) - \sin(x)\sin(y))，我们可以推导出以下恒等式
[\operatorname{versin}(x + y) = \operatorname{versin}(x)\operatorname{versin}(y) - \operatorname{versin}(x) - \operatorname{versin}(y) - \sin(x)\sin(y) + 1.]
假设我们决定在
eval()中自动展开这个：>>> from sympy import Add, sin >>> class versin(Function): ... @classmethod ... def eval(cls, x): ... # !! Not actually a good eval() method !! ... if isinstance(x, Add): ... a, b = x.as_two_terms() ... return (versin(a)*versin(b) - versin(a) - versin(b) ... - sin(a)*sin(b) + 1)这种方法递归地将
Add项分为两部分，并应用上述恒等式。>>> x, y, z = symbols('x y z') >>> versin(x + y) -sin(x)*sin(y) + versin(x)*versin(y) - versin(x) - versin(y) + 1但现在无法表示
versin(x + y)而不进行展开。这也会影响其他方法。例如，假设我们定义了微分（见下文）：>>> class versin(Function): ... @classmethod ... def eval(cls, x): ... # !! Not actually a good eval() method !! ... if isinstance(x, Add): ... a, b = x.as_two_terms() ... return (versin(a)*versin(b) - versin(a) - versin(b) ... - sin(a)*sin(b) + 1) ... ... def fdiff(self, argindex=1): ... return sin(self.args[0])我们期望
versin(x + y).diff(x)返回sin(x + y)，确实，如果我们没有在eval()中展开这个身份，它会。但使用这个版本，versin(x + y)在调用diff()之前会自动展开，因此我们得到一个更复杂的表达式：>>> versin(x + y).diff(x) sin(x)*versin(y) - sin(x) - sin(y)*cos(x)事情甚至比那更糟。让我们尝试一个有三项的
Add：>>> versin(x + y + z) (-sin(y)*sin(z) + versin(y)*versin(z) - versin(y) - versin(z) + 1)*versin(x) - sin(x)*sin(y + z) + sin(y)*sin(z) - versin(x) - versin(y)*versin(z) + versin(y) + versin(z)我们可以看到事情很快就变得失控。实际上，
versin(Add(*symbols('x:100')))（在具有 100 个项的Add上的versin()）需要超过一秒的时间来评估，而这只是创建表达式，甚至还没有进行任何操作。像这样的身份识别最好不要包含在
eval中，而是在其他方法中实现（在这种身份识别的情况下，expand_trig()）。 -
在限制输入域时：允许
None输入假设。我们的示例函数(\operatorname{versin}(x))是从(\mathbb{C})到(\mathbb{C})的函数，因此它可以接受任何输入。但假设我们有一个只对某些输入有意义的函数。作为第二个示例，让我们定义一个函数
divides如下：[\begin{split}\operatorname{divides}(m, n) = \begin{cases} 1 & \text{if}: m \mid n \ 0 & \text{if}: m\not\mid n \end{cases}.\end{split}]
也就是说，如果
m能整除n，divides(m, n)将为1，否则为0。显然，divides只在m和n为整数时有意义。我们可能会尝试像这样定义
divides的eval()方法：>>> class divides(Function): ... @classmethod ... def eval(cls, m, n): ... # !! Not actually a good eval() method !! ... ... # Evaluate for explicit integer m and n. This part is fine. ... if isinstance(m, Integer) and isinstance(n, Integer): ... return int(n % m == 0) ... ... # For symbolic arguments, require m and n to be integer. ... # If we write the logic this way, we will run into trouble. ... if not m.is_integer or not n.is_integer: ... raise TypeError("m and n should be integers")这里的问题是，通过使用
if not m.is_integer，我们要求m.is_integer必须为True。如果它是None，它将失败（有关假设为None的详细信息，请参见布尔值和三值逻辑指南）。这有两个问题。首先，它强制用户对任何输入变量定义假设。如果用户省略它们，它将失败：>>> n, m = symbols('n m') >>> print(n.is_integer) None >>> divides(m, n) Traceback (most recent call last): ... TypeError: m and n should be integers相反，他们必须编写
>>> n, m = symbols('n m', integer=True) >>> divides(m, n) divides(m, n)这似乎是一个可以接受的限制，但存在更大的问题。有时，SymPy 的假设系统无法推导出一个假设，即使在数学上是正确的。在这种情况下，它会返回
None（在 SymPy 的假设中，None表示“未定义”和“无法计算”）。例如>>> # n and m are still defined as integer=True as above >>> divides(2, (m**2 + m)/2) Traceback (most recent call last): ... TypeError: m and n should be integers在这里，表达式
(m**2 + m)/2始终是一个整数，但 SymPy 的假设系统无法推导出这一点：>>> print(((m**2 + m)/2).is_integer) NoneSymPy 的假设系统在不断改进，但总会有这样的情况，它无法推导出，这是由于问题的基本计算复杂性，以及一般问题通常是不可判定的。
因此，人们应该始终测试否定的输入变量假设，即，如果假设为
False则失败，但允许假设为None。>>> class divides(Function): ... @classmethod ... def eval(cls, m, n): ... # Evaluate for explicit integer m and n. This part is fine. ... if isinstance(m, Integer) and isinstance(n, Integer): ... return int(n % m == 0) ... ... # For symbolic arguments, require m and n to be integer. ... # This is the better way to write this logic. ... if m.is_integer is False or n.is_integer is False: ... raise TypeError("m and n should be integers")这仍然不允许非整数输入，如期望的那样：
>>> divides(1.5, 1) Traceback (most recent call last): ... TypeError: m and n should be integers但在假设为
None的情况下并不会失败：>>> divides(2, (m**2 + m)/2) divides(2, m**2/2 + m/2) >>> _.subs(m, 2) 0 >>> n, m = symbols('n m') # Redefine n and m without the integer assumption >>> divides(m, n) divides(m, n)注意
此规则仅适用于仅在引发异常时才会发生的情况，例如对输入域进行类型检查。在进行简化或其他操作的情况下，应将
None假设视为“可以是True或False”，并且不要执行可能在数学上无效的操作。 ### 假设
接下来你可能想要定义的是我们函数的假设。假设系统允许根据其输入定义函数具有的数学属性，例如，“当(x)是实数时，(f(x))是正数。”
假设系统指南详细介绍了假设系统。建议首先阅读该指南，以了解不同的假设含义以及假设系统的工作原理。
最简单的情况是一个函数始终具有给定的假设，而不考虑其输入。在这种情况下，可以直接在类上定义is_*assumption*。
例如，我们的例子divides函数总是一个整数，因为它的值总是 0 或 1：
>>> class divides(Function):
... is_integer = True
... is_negative = False
>>> divides(m, n).is_integer
True
>>> divides(m, n).is_nonnegative
True
然而，一般来说，函数的假设取决于其输入的假设。在这种情况下，应该定义一个_eval_*assumption*方法。
对于我们的(\operatorname{versin}(x))示例，当(x)是实数时，该函数始终在([0, 2])内，并且当(x)是(\pi)的偶数倍时，它恰好为 0。因此，无论x是实数还是不是π的偶数倍，versin(x)应该是非负的。记住，默认情况下，函数的定义域是(\mathbb{C})的全体，实际上versin(x)对于非实数的x也是有意义的。
要查看x是否是(\pi)的偶数倍，我们可以使用as_independent()来将x结构化地匹配为coeff*pi。在假设处理程序中，像这样结构化地分解子表达式比使用(x/pi).is_even之类的方法更可取，因为后者会创建一个新的表达式x/pi，而创建新表达式会慢得多。此外，每当创建一个表达式时，构造函数通常会导致假设被查询。如果不小心，这可能导致无限递归。因此，假设处理程序的一个好的一般规则是，永远不要在假设处理程序中创建新的表达式。始终使用像as_independent这样的结构方法来分解函数的参数。
注意(\operatorname{versin}(x))对于非实数(x)可以是非负的，例如：
>>> from sympy import I
>>> 1 - cos(pi + I*pi)
1 + cosh(pi)
>>> (1 - cos(pi + I*pi)).evalf()
12.5919532755215
对于 _eval_is_nonnegative 处理程序，如果 x.is_real 为 True，我们希望返回 True，但如果 x.is_real 为 False 或 None，则返回 None。读者可以自行处理对于使 versin(x) 非负的非实数 x 的情况，使用类似于 _eval_is_positive 处理程序的逻辑。
在假设处理方法中，就像所有方法一样，我们可以使用 self.args 访问函数的参数。
>>> from sympy.core.logic import fuzzy_and, fuzzy_not
>>> class versin(Function):
... def _eval_is_nonnegative(self):
... # versin(x) is nonnegative if x is real
... x = self.args[0]
... if x.is_real is True:
... return True
...
... def _eval_is_positive(self):
... # versin(x) is positive if x is real and not an even multiple of pi
... x = self.args[0]
...
... # x.as_independent(pi, as_Add=False) will split x as a Mul of the
... # form coeff*pi
... coeff, pi_ = x.as_independent(pi, as_Add=False)
... # If pi_ = pi, x = coeff*pi. Otherwise x is not (structurally) of
... # the form coeff*pi.
... if pi_ == pi:
... return fuzzy_and([x.is_real, fuzzy_not(coeff.is_even)])
... elif x.is_real is False:
... return False
... # else: return None. We do not know for sure whether x is an even
... # multiple of pi
>>> versin(1).is_nonnegative
True
>>> versin(2*pi).is_positive
False
>>> versin(3*pi).is_positive
True
注意在更复杂的 _eval_is_positive() 处理程序中使用 fuzzy_ 函数，并且对 if/elif 的谨慎处理很重要。在处理假设时，始终要小心正确处理三值逻辑，以确保方法在 x.is_real 或 coeff.is_even 为 None 时返回正确的答案。
警告
永远不要将 is_*assumption* 定义为 @property 方法。这样做会破坏其他假设的自动推导。is_*assumption* 应该只作为等于 True 或 False 的类变量定义。如果假设依赖于函数的 .args，则定义 _eval_*assumption* 方法。
在此示例中，不需要定义 _eval_is_real()，因为它可以从其他假设中自动推导出来，因为 nonnegative -> real。一般而言，应避免定义假设，假设系统可以根据其已知事实自动推导出的。
>>> versin(1).is_real
True
假设系统通常能够推导出比您认为的更多内容。例如，可以从上面的内容推导出当 n 是整数时，versin(2*n*pi) 为零。
>>> n = symbols('n', integer=True)
>>> versin(2*n*pi).is_zero
True
在手动编码之前，始终值得检查假设系统是否可以自动推导出某些内容。
最后，一个警告：在编写假设时，非常注意正确性。确保使用各种假设的确切定义，并始终检查是否正确处理了模糊的三值逻辑函数的 None 情况。不正确或不一致的假设可能导致微妙的错误。建议在函数具有非平凡假设处理程序时使用单元测试来检查所有不同的情况。SymPy 自身定义的所有函数都需要进行广泛测试。### 使用 evalf() 进行数值评估
这里我们展示了如何定义函数在数值上如何评估为浮点数 Float 值，例如通过 evalf()。实现数值评估可以在 SymPy 中启用多种行为。例如，一旦定义了 evalf()，您可以绘制函数，并且不等式可以评估为显式值。
如果你的函数与mpmath中的函数同名（这是 SymPy 包含的大多数函数的情况），数值评估将自动发生，你不需要做任何操作。
如果不是这种情况，可以通过定义方法_eval_evalf(self, prec)来指定数值评估，其中prec是输入的二进制精度。该方法应返回按给定精度评估的表达式，如果不可能，则返回None。
注意
_eval_evalf()方法的prec参数是二进制精度，即浮点表示中的比特数。这与evalf()方法的第一个参数不同，后者是十进制精度，即dps。例如，Float的默认二进制精度是 53，对应于十进制精度 15。因此，如果你的_eval_evalf()方法递归地调用另一个表达式的 evalf，应该调用expr._eval_evalf(prec)而不是expr.evalf(prec)，因为后者会错误地使用prec作为十进制精度。
我们可以通过递归评估(2\sin²\left(\frac{x}{2}\right))，为我们的示例(\operatorname{versin}(x))函数定义数值评估，这是编写(1 - \cos(x))更为稳定的方法。
>>> from sympy import sin
>>> class versin(Function):
... def _eval_evalf(self, prec):
... return (2*sin(self.args[0]/2)**2)._eval_evalf(prec)
>>> versin(1).evalf()
0.459697694131860
一旦定义了_eval_evalf()，就可以自动评估浮点输入。在eval()中手动实现这一点是不必要的。
>>> versin(1.)
0.459697694131860
注意evalf()可能会传递任何表达式，而不仅仅是可以数值化评估的表达式。在这种情况下，预计会对表达式的数值部分进行评估。一个常见的模式是在函数的参数上递归调用_eval_evalf(prec)。
在可能的情况下，最好重用现有 SymPy 函数中定义的 evalf 功能。但在某些情况下，需要直接使用 mpmath。 ### 重写和简化
各种简化函数和方法允许在自定义子类上指定它们的行为。并非每个 SymPy 函数都有这样的钩子。查看每个单独函数的文档以获取详细信息。
rewrite()
rewrite()方法允许根据特定函数或规则将表达式重写。例如，
>>> sin(x).rewrite(cos)
cos(x - pi/2)
要实现重写，定义一个方法_eval_rewrite(self, rule, args, **hints)，其中
-
rule是传递给rewrite()方法的规则。通常rule将是要重写为的对象的类，尽管对于更复杂的重写，它可以是任何东西。定义_eval_rewrite()的每个对象都会定义它支持的规则。许多 SymPy 函数重写为常见类，例如expr.rewrite(Add)，以执行简化或其他计算。 -
args是用于重写函数的参数。这应该使用self.args而不是self.args，因为参数中的任何递归表达式将在args中重写（假设调用者使用了rewrite(deep=True)，这是默认值）。 -
**hints是额外的关键字参数，可能用于指定重写行为。未知的提示应该被忽略，因为它们可能传递给其他_eval_rewrite()方法。如果递归调用重写，应该通过传递**hints。
该方法应返回重写的表达式，使用 args 作为函数的参数，如果表达式不应更改，则返回 None。
对于我们的 versin 示例，我们可以实现一个明显的重写，将 versin(x) 重写为 1 - cos(x)：
>>> class versin(Function):
... def _eval_rewrite(self, rule, args, **hints):
... if rule == cos:
... return 1 - cos(*args)
>>> versin(x).rewrite(cos)
1 - cos(x)
一旦我们定义了这个，simplify() 现在可以简化一些包含 versin 的表达式：
>>> from sympy import simplify
>>> simplify(versin(x) + cos(x))
1
``` #### `doit()`
`doit()` 方法 用于评估“未评估”的函数。要定义 `doit()`，实现 `doit(self, deep=True, **hints)`。如果 `deep=True`，`doit()` 应递归调用参数的 `doit()`。`**hints` 将是传递给用户的任何其他关键字参数，应该传递给 `doit()` 的任何递归调用。您可以使用 `hints` 允许用户指定 `doit()` 的特定行为。
自定义 `Function` 子类中 `doit()` 的典型用法是执行更高级的评估，这在 `eval()` 中不执行。
例如，对于我们的 `divides` 示例，有几个实例可以使用一些身份简化。例如，我们定义了 `eval()` 来评估显式整数，但我们可能也希望评估类似 `divides(k, k*n)` 这样的例子，其中除法在符号上是真实的。`eval()` 的最佳实践之一是避免过多的自动评估。在这种情况下自动评估可能被认为是过多的，因为它会使用假设系统，这可能是昂贵的。此外，我们可能希望能够表示 `divides(k, k*n)` 而不总是评估它。
解决方案是在 `doit()` 中实现这些更高级的评估。这样，我们可以通过调用 `expr.doit()` 显式执行它们，但默认情况下不会发生。例如，为 `divides` 编写的 `doit()` 可以执行这种简化（与上述的 `eval()` 定义）可能看起来像这样：
注意
如果 `doit()` 返回一个 Python `int` 文字，则将其转换为 `Integer`，以便返回的对象是 SymPy 类型。
```py
>>> from sympy import Integer
>>> class divides(Function):
... # Define evaluation on basic inputs, as well as type checking that the
... # inputs are not nonintegral.
... @classmethod
... def eval(cls, m, n):
... # Evaluate for explicit integer m and n.
... if isinstance(m, Integer) and isinstance(n, Integer):
... return int(n % m == 0)
...
... # For symbolic arguments, require m and n to be integer.
... if m.is_integer is False or n.is_integer is False:
... raise TypeError("m and n should be integers")
...
... # Define doit() as further evaluation on symbolic arguments using
... # assumptions.
... def doit(self, deep=False, **hints):
... m, n = self.args
... # Recursively call doit() on the args whenever deep=True.
... # Be sure to pass deep=True and **hints through here.
... if deep:
... m, n = m.doit(deep=deep, **hints), n.doit(deep=deep, **hints)
...
... # divides(m, n) is 1 iff n/m is an integer. Note that m and n are
... # already assumed to be integers because of the logic in eval().
... isint = (n/m).is_integer
... if isint is True:
... return Integer(1)
... elif isint is False:
... return Integer(0)
... else:
... return divides(m, n)
（注意，这使用了 约定，即 (k \mid 0) 对于所有 (k)，因此我们无需检查 m 或 n 是否为非零。如果我们使用不同的约定，我们将需要在执行简化之前检查 m.is_zero 和 n.is_zero。）
>>> n, m, k = symbols('n m k', integer=True)
>>> divides(k, k*n)
divides(k, k*n)
>>> divides(k, k*n).doit()
1
另一种常见的 doit() 实现方式是始终返回另一个表达式。这实际上将函数视为另一个表达式的“未评估”形式。
例如，我们可以定义一个 融合乘加 的函数：(\operatorname{FMA}(x, y, z) = xy + z)。将此函数表达为一个独立的函数可能对代码生成有用，但在某些情况下，将 FMA(x, y, z) “评估” 为 x*y + z 也可能很有用，以便能够与其他表达式正确简化。
>>> from sympy import Number
>>> class FMA(Function):
... """
... FMA(x, y, z) = x*y + z
... """
... @classmethod
... def eval(cls, x, y, z):
... # Number is the base class of Integer, Rational, and Float
... if all(isinstance(i, Number) for i in [x, y, z]):
... return x*y + z
...
... def doit(self, deep=True, **hints):
... x, y, z = self.args
... # Recursively call doit() on the args whenever deep=True.
... # Be sure to pass deep=True and **hints through here.
... if deep:
... x = x.doit(deep=deep, **hints)
... y = y.doit(deep=deep, **hints)
... z = z.doit(deep=deep, **hints)
... return x*y + z
>>> x, y, z = symbols('x y z')
>>> FMA(x, y, z)
FMA(x, y, z)
>>> FMA(x, y, z).doit()
x*y + z
大多数自定义函数不希望以这种方式定义 doit()。然而，这可以在始终评估的函数和从不评估的函数之间提供一个折中，从而产生一个默认情况下不评估但可以按需评估的函数（参见上文的 讨论）。 #### expand()
expand() 函数以各种方式“扩展”表达式。它实际上是几个子扩展提示的包装器。每个函数对应于 expand() 函数/方法的一个提示。可以通过定义 _eval_expand_*hint*(self, **hints) 在自定义函数中定义特定的扩展 hint。有关定义的提示以及每个特定 expand_*hint*() 函数的文档，请参阅 expand() 的文档。
**hints 关键字参数是可以传递给 expand 函数以指定额外行为的额外提示（这些与前一段描述的预定义 hints 是分开的）。未知的提示应该被忽略，因为它们可能适用于其他函数的自定义 expand() 方法。定义一个常见的提示是 force，其中 force=True 将强制进行扩展，这可能对于所有给定的输入假设在数学上并不总是有效。例如，expand_log(log(x*y), force=True) 产生 log(x) + log(y)，尽管这个恒等式并不对所有复数 x 和 y 都成立（通常 force=False 是默认值）。
注意，expand() 方法会自动处理使用其自身 deep 标志递归扩展表达式，因此 _eval_expand_* 方法不应在函数的参数上递归调用 expand。
对于我们的versin示例，我们可以通过定义一个_eval_expand_trig方法来定义trig的基本展开，该方法在1 - cos(x)上递归调用expand_trig()：
>>> from sympy import expand_trig
>>> y = symbols('y')
>>> class versin(Function):
... def _eval_expand_trig(self, **hints):
... x = self.args[0]
... return expand_trig(1 - cos(x))
>>> versin(x + y).expand(trig=True)
sin(x)*sin(y) - cos(x)*cos(y) + 1
更复杂的实现可能会尝试将expand_trig(1 - cos(x))的结果重新转换为versin函数。这留给读者作为一个练习。### 微分
要通过diff()定义微分，请定义一个方法fdiff(self, argindex)。fdiff()应该返回函数的导数，不考虑链式法则，关于第argindex个变量。argindex从1开始索引。
也就是说，f(x1, ..., xi, ..., xn).fdiff(i)应该返回(\frac{d}{d x_i} f(x_1, \ldots, x_i, \ldots, x_n))，其中(x_k)彼此独立。diff()将自动使用fdiff()的结果应用链式法则。用户代码应该使用diff()，而不是直接调用fdiff()。
注意
Function子类应该使用fdiff()来定义微分。不是Function子类的Expr的子类需要定义_eval_derivative()。不建议在Function子类上重新定义_eval_derivative()。
对于我们的(\operatorname{versin})示例函数，导数是(\sin(x))。
>>> class versin(Function):
... def fdiff(self, argindex=1):
... # argindex indexes the args, starting at 1
... return sin(self.args[0])
>>> versin(x).diff(x)
sin(x)
>>> versin(x**2).diff(x)
2*x*sin(x**2)
>>> versin(x + y).diff(x)
sin(x + y)
作为具有多个参数的函数的示例，考虑上述定义的融合乘加（FMA）示例((\operatorname{FMA}(x, y, z) = xy + z))。
我们有
[\frac{d}{dx} \operatorname{FMA}(x, y, z) = y,][\frac{d}{dy} \operatorname{FMA}(x, y, z) = x,][\frac{d}{dz} \operatorname{FMA}(x, y, z) = 1.]
因此，FMA的fdiff()方法如下所示：
>>> from sympy import Number, symbols
>>> x, y, z = symbols('x y z')
>>> class FMA(Function):
... """
... FMA(x, y, z) = x*y + z
... """
... def fdiff(self, argindex):
... # argindex indexes the args, starting at 1
... x, y, z = self.args
... if argindex == 1:
... return y
... elif argindex == 2:
... return x
... elif argindex == 3:
... return 1
>>> FMA(x, y, z).diff(x)
y
>>> FMA(x, y, z).diff(y)
x
>>> FMA(x, y, z).diff(z)
1
>>> FMA(x**2, x + 1, y).diff(x)
x**2 + 2*x*(x + 1)
要保留一个导数未求值，应该引发sympy.core.function.ArgumentIndexError(self, argindex)。如果没有定义fdiff()，这是默认行为。这里有一个在第二个参数上具有未求值导数的例子函数(f(x, y))。
>>> from sympy.core.function import ArgumentIndexError
>>> class f(Function):
... @classmethod
... def eval(cls, x, y):
... pass
...
... def fdiff(self, argindex):
... if argindex == 1:
... return 1
... raise ArgumentIndexError(self, argindex)
>>> f(x, y).diff(x)
1
>>> f(x, y).diff(y)
Derivative(f(x, y), y)
打印
您可以使用各种打印机来定义函数的打印方式，例如string printer、pretty printers和LaTeX printer，以及各种语言的代码打印机，如C和Fortran。
在大多数情况下，您不需要定义任何打印方法。默认行为是使用它们的名称打印函数。但是，在某些情况下，我们可能希望为函数定义特殊的打印方式。
例如，对于我们之前的除法示例，我们可能希望 LaTeX 打印机打印更数学化的表达式。让我们让 LaTeX 打印机表示 divides(m, n) 为 \left [ m \middle | n \right ]，看起来像是 (\left [ m \middle | n \right ])（这里 ([P]) 是Iverson 括号，如果 (P) 成立则为 (1)，否则为 (0)）。
SymPy 对象的打印方式有两种主要方法。一种是在打印机类上定义打印机。SymPy 库中的大多数类应该使用此方法，在 sympy.printing 中的各个类上定义打印机。对于用户代码，如果您定义了自定义打印机或者您有许多自定义函数需要定义打印方式，则可能更可取。参见 自定义打印机示例 了解如何以此方式定义打印机的示例。
另一种方法是在函数类上定义打印方式。要做到这一点，首先查找要为其定义打印方式的打印机上的 printmethod 属性。这是您应该为该打印机定义的方法的名称。对于 LaTeX 打印机，LatexPrinter.printmethod 是 '_latex'。打印方法总是接受一个参数 printer。应使用 printer._print 递归打印任何其他表达式，包括函数的参数。
因此，要定义我们的 divides LaTeX 打印机，我们将在类上定义如下函数 _latex(self, printer)：
>>> from sympy import latex
>>> class divides(Function):
... def _latex(self, printer):
... m, n = self.args
... _m, _n = printer._print(m), printer._print(n)
... return r'\left [ %s \middle | %s \right ]' % (_m, _n)
>>> print(latex(divides(m, n)))
\left [ m \middle | n \right ]
有关如何定义打印机方法及一些应避免的陷阱的更多详细信息，请参见 自定义打印方法示例。最重要的是，您应始终使用 printer._print() 递归打印函数的参数，包括自定义打印机内部。
其他方法
可以在自定义函数上定义几种其他方法以指定各种行为。
inverse()
inverse(self, argindex=1) 方法可以被定义为指定函数的反函数。这由 solve() 和 solveset() 使用。argindex 参数是函数的参数，从 1 开始（类似于fdiff() 方法的相同参数名称）。
inverse() 应该返回一个函数（而不是一个表达式）作为其反函数。如果反函数比单个函数更大，则可以返回一个 lambda 函数。
inverse() 应该仅对一对一的函数进行定义。换句话说，f(x).inverse() 是 f(x) 的左逆函数。在非一对一的函数上定义 inverse() 可能导致 solve() 不会给出包含该函数表达式的所有可能解。
我们的示例 versine 函数不是一对一的（因为余弦函数不是），但它的反函数 (\operatorname{arcversin}) 是。我们可以定义它如下（使用与 SymPy 中其他反三角函数相同的命名约定）：
>>> class aversin(Function):
... def inverse(self, argindex=1):
... return versin
这使得 solve() 在 aversin(x) 上工作：
>>> from sympy import solve
>>> solve(aversin(x) - y, x)
[versin(y)]
as_real_imag()
方法as_real_imag()定义如何将函数拆分为其实部和虚部。它被各种 SymPy 函数使用，这些函数分别在表达式的实部和虚部上操作。
as_real_imag(self, deep=True, **hints) 应该返回一个包含函数的实部和虚部的二元组。也就是说，expr.as_real_imag() 返回 (re(expr), im(expr))，其中 expr == re(expr) + im(expr)*I，并且 re(expr) 和 im(expr) 是实数。
如果 deep=True，它应该在其参数上递归调用 as_real_imag(deep=True, **hints)。与doit()和 eval_expand*() 方法类似，**hints 可以是任何提示，允许用户指定方法的行为。未知提示应在递归调用中传递，以防它们适用于其他 as_real_imag() 方法。
对于我们的versin 示例，我们可以递归地使用已经定义在 1 - cos(x) 上的 as_real_imag()。
>>> class versin(Function):
... def as_real_imag(self, deep=True, **hints):
... return (1 - cos(self.args[0])).as_real_imag(deep=deep, **hints)
>>> versin(x).as_real_imag()
(-cos(re(x))*cosh(im(x)) + 1, sin(re(x))*sinh(im(x)))
定义 as_real_imag() 也会自动使expand_complex()工作。
>>> versin(x).expand(complex=True)
I*sin(re(x))*sinh(im(x)) - cos(re(x))*cosh(im(x)) + 1
各种 _eval_* 方法
SymPy 中还有许多其他函数，通过自定义 _eval_* 方法可以定义这些函数的行为，类似于上述描述的方法。有关如何定义每个方法的详细信息，请参阅特定函数的文档。## 完整示例
这里是在本指南中定义的示例函数的完整示例。有关每个方法的详细信息，请参见上面的各节。
Versine
Versine（反正弦）函数定义为
[\operatorname{versin}(x) = 1 - \cos(x).]
Versine 是一个为所有复数定义的简单函数的示例。数学定义很简单，这使得在其上定义所有上述方法变得简单（在大多数情况下，我们可以重用已定义在 1 - cos(x) 上的现有 SymPy 逻辑）。
定义
>>> from sympy import Function, cos, expand_trig, Integer, pi, sin
>>> from sympy.core.logic import fuzzy_and, fuzzy_not
>>> class versin(Function):
... r"""
... The versine function.
...
... $\operatorname{versin}(x) = 1 - \cos(x) = 2\sin(x/2)².$
...
... Geometrically, given a standard right triangle with angle x in the
... unit circle, the versine of x is the positive horizontal distance from
... the right angle of the triangle to the rightmost point on the unit
... circle. It was historically used as a more numerically accurate way to
... compute 1 - cos(x), but it is rarely used today.
...
... References
... ==========
...
... .. [1] https://en.wikipedia.org/wiki/Versine
... .. [2] https://blogs.scientificamerican.com/roots-of-unity/10-secret-trig-functions-your-math-teachers-never-taught-you/
... """
... # Define evaluation on basic inputs.
... @classmethod
... def eval(cls, x):
... # If x is an explicit integer multiple of pi, x/pi will cancel and
... # be an Integer.
... n = x/pi
... if isinstance(n, Integer):
... return 1 - (-1)**n
...
... # Define numerical evaluation with evalf().
... def _eval_evalf(self, prec):
... return (2*sin(self.args[0]/2)**2)._eval_evalf(prec)
...
... # Define basic assumptions.
... def _eval_is_nonnegative(self):
... # versin(x) is nonnegative if x is real
... x = self.args[0]
... if x.is_real is True:
... return True
...
... def _eval_is_positive(self):
... # versin(x) is positive if x is real and not an even multiple of pi
... x = self.args[0]
...
... # x.as_independent(pi, as_Add=False) will split x as a Mul of the
... # form n*pi
... coeff, pi_ = x.as_independent(pi, as_Add=False)
... # If pi_ = pi, x = coeff*pi. Otherwise pi_ = 1 and x is not
... # (structurally) of the form n*pi.
... if pi_ == pi:
... return fuzzy_and([x.is_real, fuzzy_not(coeff.is_even)])
... elif x.is_real is False:
... return False
... # else: return None. We do not know for sure whether x is an even
... # multiple of pi
...
... # Define the behavior for various simplification and rewriting
... # functions.
... def _eval_rewrite(self, rule, args, **hints):
... if rule == cos:
... return 1 - cos(*args)
... elif rule == sin:
... return 2*sin(x/2)**2
...
... def _eval_expand_trig(self, **hints):
... x = self.args[0]
... return expand_trig(1 - cos(x))
...
... def as_real_imag(self, deep=True, **hints):
... # reuse _eval_rewrite(cos) defined above
... return self.rewrite(cos).as_real_imag(deep=deep, **hints)
...
... # Define differentiation.
... def fdiff(self, argindex=1):
... return sin(self.args[0])
示例
评估：
>>> x, y = symbols('x y')
>>> versin(x)
versin(x)
>>> versin(2*pi)
0
>>> versin(1.0)
0.459697694131860
假设：
>>> n = symbols('n', integer=True)
>>> versin(n).is_real
True
>>> versin((2*n + 1)*pi).is_positive
True
>>> versin(2*n*pi).is_zero
True
>>> print(versin(n*pi).is_positive)
None
>>> r = symbols('r', real=True)
>>> print(versin(r).is_positive)
None
>>> nr = symbols('nr', real=False)
>>> print(versin(nr).is_nonnegative)
None
简化：
>>> a, b = symbols('a b', real=True)
>>> from sympy import I
>>> versin(x).rewrite(cos)
1 - cos(x)
>>> versin(x).rewrite(sin)
2*sin(x/2)**2
>>> versin(2*x).expand(trig=True)
2 - 2*cos(x)**2
>>> versin(a + b*I).expand(complex=True)
I*sin(a)*sinh(b) - cos(a)*cosh(b) + 1
微分：
>>> versin(x).diff(x)
sin(x)
解决：
（一个更一般的aversin版本将定义所有上述方法）
>>> class aversin(Function):
... def inverse(self, argindex=1):
... return versin
>>> from sympy import solve
>>> solve(aversin(x**2) - y, x)
[-sqrt(versin(y)), sqrt(versin(y))]
``` ### divides
divides 是一个由以下函数定义的函数
\[\begin{split}\operatorname{divides}(m, n) = \begin{cases} 1 & \text{如果}\: m \mid n \\ 0 & \text{如果}\: m\not\mid n \end{cases},\end{split}\]
也就是说，`divides(m, n)`当`m`除以`n`时为 1，否则为 0。它仅对整数`m`和`n`定义。为了简单起见，我们使用约定\(m \mid 0\)对所有整数\(m\)成立。
`divides`是一个仅对某些输入值（整数）定义的函数的示例。`divides`还展示了如何定义自定义打印器（`_latex()`）的示例。
#### 定义
```py
>>> from sympy import Function, Integer
>>> from sympy.core.logic import fuzzy_not
>>> class divides(Function):
... r"""
... $$\operatorname{divides}(m, n) = \begin{cases} 1 & \text{for}\: m \mid n \\ 0 & \text{for}\: m\not\mid n \end{cases}.$$
...
... That is, ``divides(m, n)`` is ``1`` if ``m`` divides ``n`` and ``0``
... if ``m`` does not divide ``n`. It is undefined if ``m`` or ``n`` are
... not integers. For simplicity, the convention is used that
... ``divides(m, 0) = 1`` for all integers ``m``.
...
... References
... ==========
...
... .. [1] https://en.wikipedia.org/wiki/Divisor#Definition
... """
... # Define evaluation on basic inputs, as well as type checking that the
... # inputs are not nonintegral.
... @classmethod
... def eval(cls, m, n):
... # Evaluate for explicit integer m and n.
... if isinstance(m, Integer) and isinstance(n, Integer):
... return int(n % m == 0)
...
... # For symbolic arguments, require m and n to be integer.
... if m.is_integer is False or n.is_integer is False:
... raise TypeError("m and n should be integers")
...
... # Define basic assumptions.
...
... # divides is always either 0 or 1.
... is_integer = True
... is_negative = False
...
... # Whether divides(m, n) is 0 or 1 depends on m and n. Note that this
... # method only makes sense because we don't automatically evaluate on
... # such cases, but instead simplify these cases in doit() below.
... def _eval_is_zero(self):
... m, n = self.args
... if m.is_integer and n.is_integer:
... return fuzzy_not((n/m).is_integer)
...
... # Define doit() as further evaluation on symbolic arguments using
... # assumptions.
... def doit(self, deep=False, **hints):
... m, n = self.args
... # Recursively call doit() on the args whenever deep=True.
... # Be sure to pass deep=True and **hints through here.
... if deep:
... m, n = m.doit(deep=deep, **hints), n.doit(deep=deep, **hints)
...
... # divides(m, n) is 1 iff n/m is an integer. Note that m and n are
... # already assumed to be integers because of the logic in eval().
... isint = (n/m).is_integer
... if isint is True:
... return Integer(1)
... elif isint is False:
... return Integer(0)
... else:
... return divides(m, n)
...
... # Define LaTeX printing for use with the latex() function and the
... # Jupyter notebook.
... def _latex(self, printer):
... m, n = self.args
... _m, _n = printer._print(m), printer._print(n)
... return r'\left [ %s \middle | %s \right ]' % (_m, _n)
...
示例
评估
>>> from sympy import symbols
>>> n, m, k = symbols('n m k', integer=True)
>>> divides(3, 10)
0
>>> divides(3, 12)
1
>>> divides(m, n).is_integer
True
>>> divides(k, 2*k)
divides(k, 2*k)
>>> divides(k, 2*k).is_zero
False
>>> divides(k, 2*k).doit()
1
打印：
>>> str(divides(m, n)) # This is using the default str printer
'divides(m, n)'
>>> print(latex(divides(m, n)))
\left [ m \middle | n \right ]
``` ### 融合乘加（FMA）
[融合乘加（FMA）](https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation#Fused_multiply%E2%80%93add)是一种先乘后加：
\[\operatorname{FMA}(x, y, z) = xy + z.\]
它通常在硬件上实现为单一浮点操作，具有比乘法和加法组合更好的舍入和性能。
FMA 是一个自定义函数的示例，它被定义为另一个函数的未评估的“缩写”。这是因为`doit()`方法被定义为返回`x*y + z`，这意味着`FMA`函数可以轻松地评估为它代表的表达式，但`eval()`方法并不返回任何内容（除非`x`、`y`和`z`都是明确的数值），这意味着默认情况下它保持未评估状态。
与 versine 示例相比，它将`versin`视为自己的一级函数。尽管`versin(x)`可以用其他函数（`1 - cos(x)`）来表达，但在`versin.eval()`中不会对一般的符号输入进行评估，而且`versin.doit()`根本没有定义。
`FMA`也是一个在多个变量上定义的连续函数的示例，它展示了在`fdiff`示例中`argindex`的工作方式。
最后，`FMA`展示了为`C`和`C++`定义一些代码打印器的示例（使用来自`C99CodePrinter.printmethod`和`CXX11CodePrinter.printmethod`的方法名称），因为这是该函数的典型用例。
FMA 的数学定义非常简单，定义每种方法都很容易，但这里只展示了少数几种。正弦和除法示例展示了如何定义本指南讨论的其他重要方法。
请注意，如果您想要实际在代码生成中使用融合乘加法，SymPy 中已经有一个版本`sympy.codegen.cfunctions.fma()`，它受现有代码打印机的支持。这里的版本仅设计为示例。
#### 定义
```py
>>> from sympy import Number, symbols, Add, Mul
>>> x, y, z = symbols('x y z')
>>> class FMA(Function):
... """
... FMA(x, y, z) = x*y + z
...
... FMA is often defined as a single operation in hardware for better
... rounding and performance.
...
... FMA can be evaluated by using the doit() method.
...
... References
... ==========
...
... .. [1] https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation#Fused_multiply%E2%80%93add
... """
... # Define automatic evaluation on explicit numbers
... @classmethod
... def eval(cls, x, y, z):
... # Number is the base class of Integer, Rational, and Float
... if all(isinstance(i, Number) for i in [x, y, z]):
... return x*y + z
...
... # Define numerical evaluation with evalf().
... def _eval_evalf(self, prec):
... return self.doit(deep=False)._eval_evalf(prec)
...
... # Define full evaluation to Add and Mul in doit(). This effectively
... # treats FMA(x, y, z) as just a shorthand for x*y + z that is useful
... # to have as a separate expression in some contexts and which can be
... # evaluated to its expanded form in other contexts.
... def doit(self, deep=True, **hints):
... x, y, z = self.args
... # Recursively call doit() on the args whenever deep=True.
... # Be sure to pass deep=True and **hints through here.
... if deep:
... x = x.doit(deep=deep, **hints)
... y = y.doit(deep=deep, **hints)
... z = z.doit(deep=deep, **hints)
... return x*y + z
...
... # Define FMA.rewrite(Add) and FMA.rewrite(Mul).
... def _eval_rewrite(self, rule, args, **hints):
... x, y, z = self.args
... if rule in [Add, Mul]:
... return self.doit()
...
... # Define differentiation.
... def fdiff(self, argindex):
... # argindex indexes the args, starting at 1
... x, y, z = self.args
... if argindex == 1:
... return y
... elif argindex == 2:
... return x
... elif argindex == 3:
... return 1
...
... # Define code printers for ccode() and cxxcode()
... def _ccode(self, printer):
... x, y, z = self.args
... _x, _y, _z = printer._print(x), printer._print(y), printer._print(z)
... return "fma(%s, %s, %s)" % (_x, _y, _z)
...
... def _cxxcode(self, printer):
... x, y, z = self.args
... _x, _y, _z = printer._print(x), printer._print(y), printer._print(z)
... return "std::fma(%s, %s, %s)" % (_x, _y, _z)
示例
评估：
>>> x, y, z = symbols('x y z')
>>> FMA(2, 3, 4)
10
>>> FMA(x, y, z)
FMA(x, y, z)
>>> FMA(x, y, z).doit()
x*y + z
>>> FMA(x, y, z).rewrite(Add)
x*y + z
>>> FMA(2, pi, 1).evalf()
7.28318530717959
微分
>>> FMA(x, x, y).diff(x)
2*x
>>> FMA(x, y, x).diff(x)
y + 1
代码打印机
>>> from sympy import ccode, cxxcode
>>> ccode(FMA(x, y, z))
'fma(x, y, z)'
>>> cxxcode(FMA(x, y, z))
'std::fma(x, y, z)'
附加提示
-
SymPy 包含数十个函数。这些可以作为编写自定义函数的有用示例，特别是如果函数类似于已实现的函数。请记住，本指南中的所有内容同样适用于 SymPy 中包含的函数和用户定义的函数。实际上，本指南旨在既是 SymPy 贡献者的开发指南，也是 SymPy 最终用户的指南。
-
如果您有许多共享常规逻辑的自定义函数，可以使用一个通用的基类来包含这些共享逻辑。例如，请参阅SymPy 中三角函数的源代码，其中使用了
TrigonometricFunction、InverseTrigonometricFunction和ReciprocalTrigonometricFunction基类及其一些共享逻辑。 -
与任何代码一样，为您的函数编写广泛的测试是一个好主意。SymPy 测试套件提供了有关如何为这些函数编写测试的示例。SymPy 本身包含的所有代码都必须进行测试。SymPy 包含的函数还应始终包含带有引用、数学定义和 doctest 示例的文档字符串。
物理学
Python 包 SymPy 可用于符号解决经典力学、连续介质力学、控制系统、高能物理、光学、量子力学、单位和向量代数等物理问题。
- 控制包示例
控制包示例
原文：
docs.sympy.org/latest/guides/physics/control_problems.html
下面提供了一些全面的教材示例，演示控制模块可能的用例。
示例 1
给出未知传递函数的极零图如上所示。
-
确定系统的连续时间DC 增益为20时的确切传递函数。
-
传递函数稳定还是不稳定的性质。
-
获取系统的单位冲击响应。
-
找到系统的时域响应的初始值，不使用时域方程。
解决方案
>>> # Imports >>> from sympy import symbols, I, limit, pprint, solve, oo >>> from sympy.physics.control import TransferFunction子部分 1
>>> s, k = symbols('s k') >>> gain = k # Let unknwon gain be k >>> a = [-3] # Zero at -3 in S plane >>> b = [-1, -2-I, -2+I] # Poles at -1, (-2, j) and (-2, -j) in S plane >>> tf = TransferFunction.from_zpk(a, b, gain, s) >>> pprint(tf) k*(s + 3) ------------------------------- (s + 1)*(s + 2 - I)*(s + 2 + I) >>> gain = tf.dc_gain() >>> print(gain) 3*k*(2 - I)*(2 + I)/25 >>> K = solve(gain - 20, k)[0] # Solve for k >>> tf = tf.subs({k: K}) # Reconstruct the TransferFunction using .subs() >>> pprint(tf.expand()) 100*s ----- + 100 3 ------------------- 3 2 s + 5*s + 9*s + 5子部分 2
>>> tf.is_stable() # Expect True, since poles lie in the left half of S plane True子部分 3
>>> from sympy import inverse_laplace_transform >>> t = symbols('t', positive = True) >>> # Convert from S to T domain for impulse response >>> tf = tf.to_expr() >>> Impulse_Response = inverse_laplace_transform(tf, s, t) >>> pprint(Impulse_Response) -t -2*t 100*e 100*e *cos(t) ------- - ---------------- 3 3子部分 4
>>> # Apply the Initial Value Theorem on Equation of S domain >>> # limit(y(t), t, 0) = limit(s*Y(S), s, oo) >>> limit(s*tf, s, oo) 0
示例 2
找到以下弹簧-质量-阻尼系统的传递函数：
解决方案
>>> # Imports
>>> from sympy import Function, laplace_transform, laplace_initial_conds, laplace_correspondence, diff, Symbol, solve
>>> from sympy.abc import s, t
>>> from sympy.physics.control import TransferFunction
>>> y = Function('y')
>>> Y = Function('Y')
>>> u = Function('u')
>>> U = Function('U')
>>> k = Symbol('k') # Spring Constant
>>> c = Symbol('c') # Damper
>>> m = Symbol('m') # Mass of block
系统的微分方程如下所示：
[\begin{split}\frac{{d²y(t)}}{{dt²}} + c\frac{{dy(t)}}{{dt}} + ky(t) = w²u(t) \\ with \ initial \ conditions \ y(0) = t,\quad\frac{{dy}}{{dt}}\bigg|_{t=0} = 0\\end{split}]
>>> f = m*diff(y(t), t, t) + c*diff(y(t), t) + k*y(t) - u(t) >>> F = laplace_transform(f, t, s, noconds=True) >>> F = laplace_correspondence(F, {u: U, y: Y}) >>> F = laplace_initial_conds(F, t, {y: [0, 0]}) >>> t = (solve(F, Y(s))[0])/U(s) # To construct Transfer Function from Y(s) and U(s) >>> tf = TransferFunction.from_rational_expression(t, s) >>> pprint(tf) 1 -------------- 2 c*s + k + m*s
示例 3
在时间域中称为冲激响应矩阵的信号矩阵 g(t) 如下所示。
[\begin{split}g(t) = \begin{bmatrix} (1-t)e^{-t} & e^{-2t} \ -e{-t}+5e & \left(-3\sqrt{3}\sin\left(\frac{\sqrt{3}t}{2}\right)+\cos\left(\frac{\sqrt{3}t}{2}\right)\right)e^{-\frac{t}{2}} \end{bmatrix}\end{split}]
关于此矩阵，找到
-
系统矩阵（拉普拉斯域的传递函数矩阵）(g(t) → G(s))。
-
系统中输入和输出信号的数量。
-
系统元素（传递函数矩阵中各个传递函数的极点和零点）在拉普拉斯域的极点和零点 (注意：多输入多输出系统的实际极点和零点并非传递函数矩阵中各个元素的极点和零点)。同时，可视化与G(s)矩阵的第 1 个输入和第 1 个输出对应的单个传递函数的极点和零点。
-
绘制与G(s)矩阵的第 1 个输入和第 1 个输出对应的单个传递函数的单位阶跃响应。
-
分析与G(s)矩阵的第 1 个输入和第 2 个输出对应的传递函数的博德幅度和相位图。
解决方案
>>> # Imports >>> from sympy import Matrix, laplace_transform, inverse_laplace_transform, exp, cos, sqrt, sin, pprint >>> from sympy.abc import s, t >>> from sympy.physics.control import *子部分 1
>>> g = Matrix([[exp(-t)*(1 - t), exp(-2*t)], [5*exp((-2*t))-exp((-t)), (cos((sqrt(3)*t)/2) - 3*sqrt(3)*sin((sqrt(3)*t)/2))*exp(-t/2)]]) >>> G = g.applyfunc(lambda a: laplace_transform(a, t, s)[0]) >>> pprint(G) [ 1 1 1 ] [----- - -------- ----- ] [s + 1 2 s + 2 ] [ (s + 1) ] [ ] [ 5 1 s + 1/2 9 ] [ ----- - ----- -------------- - ------------------] [ s + 2 s + 1 2 3 / 2 3\] [ (s + 1/2) + - 2*|(s + 1/2) + -|] [ 4 \ 4/]子部分 2
>>> G = TransferFunctionMatrix.from_Matrix(G, s) >>> type(G) <class 'sympy.physics.control.lti.TransferFunctionMatrix'> >>> type(G[0]) <class 'sympy.physics.control.lti.TransferFunction'> >>> print(f'Inputs = {G.num_inputs}, Outputs = {G.num_outputs}') Inputs = 2, Outputs = 2子部分 3
>>> G.elem_poles() [[[-1, -1, -1], [-2]], [[-2, -1], [-1/2 - sqrt(3)*I/2, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]]] >>> G.elem_zeros() [[[-1, 0], []], [[-3/4], [4, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]]] >>> pole_zero_plot(G[0, 0])(
png,hires.png,
子部分 4
>>> tf1 = G[0, 0] >>> pprint(tf1) 2 -s + (s + 1) - 1 ----------------- 3 (s + 1) >>> step_response_plot(tf1)(
png,hires.png,
子部分 5
>>> tf2 = G[0, 1] >>> bode_magnitude_plot(tf2)(
png,hires.png,
>>> bode_phase_plot(tf2)(
png,hires.png,
示例 4
-
通过将P(s)和C(s)按串联配置设计一个系统，计算等效系统矩阵，当块的顺序颠倒时（即先 C(s)，然后 P(s)）。
[\begin{split}P(s) = \begin{bmatrix} \frac{1}{s} & \frac{2}{s+2} \ 0 & 3 \end{bmatrix}\end{split}][\begin{split}C(s) = \begin{bmatrix} 1 & 1 \ 2 & 2 \end{bmatrix}\end{split}]
-
同样，为了系统（负反馈环）找到等效闭环系统（或从下面给定的框图中找到 v/u 的比率），其中C(s)作为控制器，P(s)作为装置（参考下面的框图）。
![../../_images/Control_Problems_Q4.svg]()
解决方案
>>> # Imports >>> from sympy import Matrix, pprint >>> from sympy.abc import s, t >>> from sympy.physics.control import *子部分 1
>>> P_mat = Matrix([[1/s, 2/(2+s)], [0, 3]]) >>> C_mat = Matrix([[1, 1], [2, 2]]) >>> P = TransferFunctionMatrix.from_Matrix(P_mat, var=s) >>> C = TransferFunctionMatrix.from_Matrix(C_mat, var=s) >>> # Series equivalent, considering (Input)→[P]→[C]→(Output). Note that order of matrix multiplication is opposite to the order in which the elements are arranged. >>> pprint(C*P) [1 1] [1 2 ] [- -] [- -----] [1 1] [s s + 2] [ ] *[ ] [2 2] [0 3 ] [- -] [- - ] [1 1]{t} [1 1 ]{t} >>> # Series equivalent, considering (Input)→[C]→[P]→(Output). >>> pprint(P*C) [1 2 ] [1 1] [- -----] [- -] [s s + 2] [1 1] [ ] *[ ] [0 3 ] [2 2] [- - ] [- -] [1 1 ]{t} [1 1]{t} >>> pprint((C*P).doit()) [1 3*s + 8 ] [- ------- ] [s s + 2 ] [ ] [2 6*s + 16] [- --------] [s s + 2 ]{t} >>> pprint((P*C).doit()) [ 5*s + 2 5*s + 2 ] [--------- ---------] [s*(s + 2) s*(s + 2)] [ ] [ 6 6 ] [ - - ] [ 1 1 ]{t}子部分 2
>>> tfm_feedback = MIMOFeedback(P, C, sign=-1) >>> pprint(tfm_feedback.doit()) # ((I + P*C)**-1)*P [ 7*s + 14 -s - 6 ] [--------------- ---------------] [ 2 2 ] [7*s + 19*s + 2 7*s + 19*s + 2] [ ] [ 2 ] [ -6*s - 12 3*s + 9*s + 6 ] [--------------- ---------------] [ 2 2 ] [7*s + 19*s + 2 7*s + 19*s + 2]{t}
示例 5
给定，
[ \begin{align}\begin{aligned}\begin{split}G1 &= \frac{1}{10 + s}\\\end{split}\\begin{split}G2 &= \frac{1}{1 + s}\\\end{split}\\begin{split}G3 &= \frac{1 + s²}{4 + 4s + s²}\\\end{split}\\begin{split}G4 &= \frac{1 + s}{6 + s}\\\end{split}\\begin{split}H1 &= \frac{1 + s}{2 + s}\\\end{split}\\begin{split}H2 &= \frac{2 \cdot (6 + s)}{1 + s}\\\end{split}\\begin{split}H3 &= 1\\end{split}\end{aligned}\end{align} ]
其中(s)是传递函数（在拉普拉斯域）的变量。
找到
-
表示上述系统的等效传递函数。
-
极点零点图的系统。
解决方案
>>> from sympy.abc import s >>> from sympy.physics.control import * >>> G1 = TransferFunction(1, 10 + s, s) >>> G2 = TransferFunction(1, 1 + s, s) >>> G3 = TransferFunction(1 + s**2, 4 + 4*s + s**2, s) >>> G4 = TransferFunction(1 + s, 6 + s, s) >>> H1 = TransferFunction(1 + s, 2 + s, s) >>> H2 = TransferFunction(2*(6 + s), 1 + s, s) >>> H3 = TransferFunction(1, 1, s) >>> sys1 = Series(G3, G4) >>> sys2 = Feedback(sys1, H1, 1).doit() >>> sys3 = Series(G2, sys2) >>> sys4 = Feedback(sys3, H2).doit() >>> sys5 = Series(G1, sys4) >>> sys6 = Feedback(sys5, H3) >>> sys6 # Final unevaluated Feedback object Feedback(Series(TransferFunction(1, s + 10, s), TransferFunction((s + 1)**3*(s + 2)*(s + 6)**2*(s**2 + 1)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4)**2, (s + 1)*(s + 6)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*((s + 1)**2*(s + 6)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4) + (s + 1)*(s + 2)*(s + 6)*(2*s + 12)*(s**2 + 1)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4), s)), TransferFunction(1, 1, s), -1) >>> sys6.doit() # Reducing to TransferFunction form without simplification TransferFunction((s + 1)**4*(s + 2)*(s + 6)**3*(s + 10)*(s**2 + 1)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))**2*((s + 1)**2*(s + 6)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4) + (s + 1)*(s + 2)*(s + 6)*(2*s + 12)*(s**2 + 1)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4)**3, (s + 1)*(s + 6)*(s + 10)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*((s + 1)**2*(s + 6)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4) + (s + 1)*(s + 2)*(s + 6)*(2*s + 12)*(s**2 + 1)*(s**2 + 4*s + 4))*((s + 1)**3*(s + 2)*(s + 6)**2*(s**2 + 1)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4)**2 + (s + 1)*(s + 6)*(s + 10)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*((s + 1)**2*(s + 6)*(-(s + 1)**2*(s**2 + 1) + (s + 2)*(s + 6)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4) + (s + 1)*(s + 2)*(s + 6)*(2*s + 12)*(s**2 + 1)*(s**2 + 4*s + 4))*(s**2 + 4*s + 4))*(s**2 + 4*s + 4), s) >>> sys6 = sys6.doit(cancel=True, expand=True) # Simplified TransferFunction form >>> sys6 TransferFunction(s**4 + 3*s**3 + 3*s**2 + 3*s + 2, 12*s**5 + 193*s**4 + 873*s**3 + 1644*s**2 + 1484*s + 712, s) >>> pole_zero_plot(sys6)(
png,hires.png,
参考资料
解方程
Python 包 SymPy 可以符号性地解决方程、微分方程、线性方程、非线性方程、矩阵问题、不等式、丢番图方程和评估积分。SymPy 也可以进行数值解析。
解决指南页面提供适用于许多类型解决任务的建议。
学习如何使用 SymPy 计算代数系统来：
| 描述 | 示例 | 解决方案 |
|---|---|---|
| 代数方法解方程 | (x² = y) | (x \in {-\sqrt{y},\sqrt{y}}) |
| 代数方法解方程组 | (x² + y = 2z, y = -4z) | ({(x = -\sqrt{6z}, y = -4z),) ({(x = \sqrt{6z}, y = -4z)}}) |
| 数值方法求解方程（或方程组） | (\cos(x) = x ) | ( x \approx 0.739085133215161) |
| 代数方法求解常微分方程 | (y''(x) + 9y(x)=0 ) | ( y(x)=C_{1} \sin(3x)+ C_{2} \cos(3x)) |
| 代数方法求多项式的根（代数或数值方法） | ( ax² + bx + c = 0 ) | ( x = \frac{-b\pm\sqrt{b² - 4ac}}{2a} ) |
| 代数方法求解矩阵方程 | ( \left[\begin{array}{cc} c & d\1 & -e\end{array}\right] \left[\begin{array}{cc} x\y\end{array}\right] = \left[\begin{array}{cc} 2\0\end{array}\right] ) | ( \left[\begin{array}{cc} x\y\end{array}\right] = \left[\begin{array}{cc} \frac{2e}{ce+d}\\frac{2}{ce+d}\end{array}\right]) |
| 代数方法简化单变量不等式或不等式系统 | ( x² < \pi, x > 0 ) | ( 0 < x < \sqrt{\pi} ) |
| 代数方法求解丢番图方程 | (a² + b² = c²) | ((a=2pq, b=p²-q², c=p²+q²)) |
注释：
-
SymPy 有一个名为
solve()的函数，用于找到方程或方程组的解，或者函数的根。SymPy 的solve()可能或可能不适合您的特定问题，因此我们建议您使用本页上的链接来学习如何“解决”您的问题。 -
尽管一个常见的口头表达是例如“解决一个积分，”在 SymPy 的术语中，它将是“评估一个积分”。此页面不提供此类任务的指导。请搜索文档以找到您想要评估的表达类型。
解决指南
原文：
docs.sympy.org/latest/guides/solving/solving-guidance.html
这些准则适用于许多类型的解决方案。
数值解
没有封闭形式解的方程
绝大多数任意非线性方程都没有封闭形式解。可解类方程基本上是：
-
线性方程
-
多项式，除非受到Abel-Ruffini theorem的限制（了解使用
GroebnerBasis解决多项式的更多信息） -
可以通过反转某些超越函数来解决的方程
-
可以转换为上述情况的问题（例如，通过将三角函数转换为多项式）
-
还有一些特殊情况，可以用类似
Lambert W function解决 -
您可以通过任何上述方法
decompose()解决的方程
SymPy 可能会反映您的方程无法用代数（符号）方式表达的解，或者 SymPy 缺乏找到已存在的封闭形式解的算法，例如通过返回诸如NotImplementedError的错误：
>>> from sympy import solve, cos
>>> from sympy.abc import x
>>> solve(cos(x) - x, x, dict=True)
Traceback (most recent call last):
...
NotImplementedError: multiple generators [x, cos(x)]
No algorithms are implemented to solve equation -x + cos(x)
因此，您可能需要使用nsolve()之类的方法进行数值解决。
>>> from sympy import nsolve, cos
>>> from sympy.abc import x
>>> nsolve(cos(x) - x, x, 2)
0.739085133215161
如果您收到像CRootOf()这样的非封闭形式解（表示多项式的索引复数根），您可以使用evalf()进行数值评估：
>>> from sympy import solve
>>> from sympy.abc import x
>>> solutions = solve(x**5 - x - 1, x, dict=True)
>>> solutions
[{x: CRootOf(x**5 - x - 1, 0)}, {x: CRootOf(x**5 - x - 1, 1)}, {x: CRootOf(x**5 - x - 1, 2)}, {x: CRootOf(x**5 - x - 1, 3)}, {x: CRootOf(x**5 - x - 1, 4)}]
>>> [solution[x].evalf(3) for solution in solutions]
[1.17, -0.765 - 0.352*I, -0.765 + 0.352*I, 0.181 - 1.08*I, 0.181 + 1.08*I]
您可能更喜欢数值解的情况
即使您的问题有封闭形式解，您可能更喜欢数值解。
像solve()和solveset()这样的解函数将不会尝试找到数值解，只会找到数学上精确的符号解。因此，如果您需要数值解，考虑使用nsolve()。
在某些情况下，即使有闭合形式的解，也可能太繁琐而不可取。在这种情况下，如果接受数值解，则可以使用evalf()。例如，以下解集在精确表示时总共包含超过 40 项（如果需要查看所有内容，请在下面的代码块中水平滚动），而数值表示时只有八项：
>>> from sympy import symbols, solve
>>> x = symbols('x')
>>> solutions = solve(x**4 + 10*x**2 + x + 1, x, dict=True)
>>> solutions
[{x: -sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3))/2 - sqrt(-40/3 - 2*(1307/432 + sqrt(434607)*I/144)**(1/3) + 2/sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3)) - 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)))/2}, {x: sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3))/2 - sqrt(-40/3 - 2*(1307/432 + sqrt(434607)*I/144)**(1/3) - 2/sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3)) - 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)))/2}, {x: sqrt(-40/3 - 2*(1307/432 + sqrt(434607)*I/144)**(1/3) - 2/sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3)) - 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)))/2 + sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3))/2}, {x: sqrt(-40/3 - 2*(1307/432 + sqrt(434607)*I/144)**(1/3) + 2/sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3)) - 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)))/2 - sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3))/2}]
>>> for solution in solutions:
... solution[x].evalf()
-0.0509758447494279 + 0.313552108895239*I
0.0509758447494279 + 3.14751999969868*I
0.0509758447494279 - 3.14751999969868*I
-0.0509758447494279 - 0.313552108895239*I
在其他情况下，即使确切解具有少量项，您可能希望获得数值解，以便知道其近似数值。例如，估计(\sqrt{2} e^{\pi}/2)约为(16)可能会很困难：
>>> from sympy import pi, sqrt, exp, solve, evalf
>>> shorter = solve(sqrt(2)*x - exp(pi), x, dict=True)
>>> shorter
[{x: sqrt(2)*exp(pi)/2}]
>>> [solution[x].evalf(3) for solution in shorter]
[16.4]
使用精确值
如果你想保留诸如超越数和平方根等符号的精确数学值，请定义它们以便 SymPy 可以进行符号解释，例如使用 SymPy 的Pi：
>>> from sympy import symbols, solve, pi
>>> x = symbols('x')
>>> solve(x**2 - pi, x, dict=True)
[{x: -sqrt(pi)}, {x: sqrt(pi)}]
如果使用标准的 Python math 版本的(\pi)，Python 将传递该不精确值给 SymPy，导致一个不精确的数值解：
>>> from sympy import symbols, solve
>>> from math import pi
>>> x = symbols('x')
>>> solve(x**2 - pi, x, dict=True)
[{x: -1.77245385090552}, {x: 1.77245385090552}]
要使用像(6.2)或(1/2)这样的精确数值，请参阅 Python numbers vs. SymPy Numbers。
在某些情况下，使用不精确值将阻止 SymPy 找到结果。例如，可以解决这个精确方程：
>>> from sympy import symbols, solve, sqrt
>>> x = symbols('x')
>>> eq = x**sqrt(2) - 2
>>> solve(eq, x, dict=True)
[{x: 2**(sqrt(2)/2)}]
但如果使用不精确方程 eq = x**1.4142135623730951 - 2，尽管尝试了很长时间，SymPy 也不会返回结果。
在函数调用中包括要解的变量
我们建议您在包括solve()和solveset()等解决函数的第二个参数中包括要解的变量。虽然这对于一元方程是可选的，但这是一个良好的实践，因为它确保 SymPy 将解决所需的符号。例如，您可能对(x)的解决方案感兴趣，但 SymPy 却解决了(y)：
>>> from sympy.abc import x, y
>>> from sympy import solve
>>> solve(x**2 - y, dict=True)
[{y: x**2}]
指定要解的变量确保 SymPy 对其进行求解：
>>> from sympy.abc import x, y
>>> from sympy import solve
>>> solve(x**2 - y, x, dict=True)
[{x: -sqrt(y)}, {x: sqrt(y)}]
确保从solve()保持一致的格式化
solve()根据解的类型输出产生各种输出。使用dict=True将提供一致的输出格式，在以编程方式提取解决方案信息时尤其重要。
要提取解决方案，可以遍历字典列表：
>>> from sympy import parse_expr, solve, solveset
>>> from sympy.abc import x
>>> expr = "x² = y"
>>> parsed = parse_expr(expr, transformations="all")
>>> parsed
Eq(x**2, y)
>>> solutions = solve(parsed, x, dict=True)
>>> [solution[x] for solution in solutions]
[-sqrt(y), sqrt(y)]
>>> solveset(parsed, x)
{-sqrt(y), sqrt(y)}
``` ## 可加速 `solve()` 的选项
### 包括使任何分母为零的解
通常情况下，`solve()` 检查是否有任何解使任何分母为零，并自动排除它们。如果您希望包括这些解，并加速 `solve()`（尽管可能获得无效解），请设置 `check=False`：
```py
>>> from sympy import Symbol, sin, solve
>>> x = Symbol("x")
>>> solve(sin(x)/x, x, dict=True) # 0 is excluded
[{x: pi}]
>>> solve(sin(x)/x, x, dict=True, check=False) # 0 is not excluded
[{x: 0}, {x: pi}]
不要简化解决方案
通常情况下，solve() 在返回许多结果之前简化它们，并且（如果 check 不为 False）在解决方案和将它们代入函数应为零的表达式时使用一般的 simplify() 函数。如果您不希望简化解决方案，并希望加速 solve()，请使用 simplify=False。
>>> from sympy import solve
>>> from sympy.abc import x, y
>>> expr = x**2 - (y**5 - 3*y**3 + y**2 - 3)
>>> solve(expr, x, dict=True)
[{x: -sqrt(y**5 - 3*y**3 + y**2 - 3)}, {x: sqrt(y**5 - 3*y**3 + y**2 - 3)}]
>>> solve(expr, x, dict=True, simplify=False)
[{x: -sqrt((y + 1)*(y**2 - 3)*(y**2 - y + 1))}, {x: sqrt((y + 1)*(y**2 - 3)*(y**2 - y + 1))}]
解析表示方程的字符串
如果您正在创建表达式本身，则建议不要使用字符串解析来创建表达式。但是，如果您以编程方式读取字符串，则此方法很方便。
您可以解析表示方程的字符串为 SymPy 可以理解的形式（例如 Eq 形式），然后解决解析后的表达式。从字符串解析方程式需要您使用 SymPy 的 transformations。
-
解释等号
-
从您的变量创建符号
-
使用更多数学（而不是标准的 Python）符号，例如指数运算符可以从
^解析，而不必使用 Python 的**。
如果您已经在Eq（等式）形式中有方程式，则可以解析该字符串：
>>> from sympy import parse_expr, solve, solveset
>>> from sympy.abc import x
>>> expr = "Eq(x², y)"
>>> parsed = parse_expr(expr, transformations="all")
>>> parsed
Eq(x**2, y)
SymPy 还可以使用 parse_latex() 解析 LaTeX 表达式。
报告 Bug
如果您发现这些命令有 bug，请在SymPy 邮件列表上发布问题。
代数方式解方程
原文：
docs.sympy.org/latest/guides/solving/solve-equation-algebraically.html
使用 SymPy 以代数方式（符号方式）解方程。例如，解决 (x² = y) 对 (x) 的方程得出 (x \in {-\sqrt{y},\sqrt{y}})。
考虑的其他选择
-
SymPy 还可以解决许多其他类型的问题，包括方程组。
-
有些方程无法代数方式解决（无论是完全还是通过 SymPy），因此您可能需要通过数值方法解方程。
解决函数
有两个高级函数用于解方程，solve()和solveset()。以下是每个的一个示例：
solve()
>>> from sympy.abc import x, y
>>> from sympy import solve
>>> solve(x**2 - y, x, dict=True)
[{x: -sqrt(y)}, {x: sqrt(y)}]
solveset()
>>> from sympy import solveset
>>> from sympy.abc import x, y
>>> solveset(x**2 - y, x)
{-sqrt(y), sqrt(y)}
以下是何时使用的建议：
-
solve()-
您希望得到变量满足方程的不同值的显式符号表示。
-
您希望将这些显式解值替换为涉及相同变量的其他方程或表达式，使用
subs()。
-
-
solveset()-
您希望以数学上准确的方式表示解，使用数学集合。
-
您希望得到所有解的表示，包括如果存在无限多解时。
-
您希望一个一致的输入接口。
-
您希望限制解的定义域为任意集。
-
您不需要从解集中以程序方式提取解：解集不能以程序方式查询。
-
指南
请参考在函数调用中包含要解决的变量和确保从 solve()获得一致的格式。
代数方式解方程
你可以用几种方式解方程。以下示例演示了在适用的情况下同时使用 solve() 和 solveset()。你可以选择最适合你方程的函数。
将你的方程转化为等于零的表达式。
使用这样一个事实：任何不在 Eq（等式）中的表达式都会被解函数自动假定为等于零（0）。你可以将方程 (x² = y) 重新排列为 (x² - y = 0)，然后解决这个表达式。如果你正在交互地解决一个已经等于零的表达式，或者一个你不介意重新排列成 (expression = 0) 的方程，这种方法就很方便。
>>> from sympy import solve, solveset
>>> from sympy.abc import x, y
>>> solve(x**2 - y, x, dict=True)
[{x: -sqrt(y)}, {x: sqrt(y)}]
>>> solveset(x**2 - y, x)
{-sqrt(y), sqrt(y)}
将你的方程放入 Eq 形式中。
将你的方程放入 Eq 形式中，然后解 Eq。如果你正在交互地解决一个你已经有了等式形式的方程，或者你将其视为等式的方程，这种方法很方便。它还有助于在从一边减去另一边时避免符号错误。
>>> from sympy import Eq, solve, solveset
>>> from sympy.abc import x, y
>>> eqn = Eq(x**2, y)
>>> eqn
Eq(x**2, y)
>>> solutions = solve(eqn, x, dict=True)
>>> print(solutions)
[{x: -sqrt(y)}, {x: sqrt(y)}]
>>> solutions_set = solveset(eqn, x)
>>> print(solutions_set)
{-sqrt(y), sqrt(y)}
>>> for solution_set in solutions_set:
... print(solution_set)
sqrt(y)
-sqrt(y)
限制解的域。
默认情况下，SymPy 将在复数域中返回解，这也包括纯实数和纯虚数值。这里，前两个解是实数，最后两个是虚数：
>>> from sympy import Symbol, solve, solveset
>>> x = Symbol('x')
>>> solve(x**4 - 256, x, dict=True)
[{x: -4}, {x: 4}, {x: -4*I}, {x: 4*I}]
>>> solveset(x**4 - 256, x)
{-4, 4, -4*I, 4*I}
要将返回的解限制为实数，或者另一个域或范围，不同的解函数使用不同的方法。
对于 solve()，在要解的符号 (x) 上放置一个假设，
>>> from sympy import Symbol, solve
>>> x = Symbol('x', real=True)
>>> solve(x**4 - 256, x, dict=True)
[{x: -4}, {x: 4}]
或者使用标准的 Python 过滤列表技术来限制解，例如列表推导式：
>>> from sympy import Or, Symbol, solve
>>> x = Symbol('x', real=True)
>>> expr = (x-4)*(x-3)*(x-2)*(x-1)
>>> solution = solve(expr, x)
>>> print(solution)
[1, 2, 3, 4]
>>> solution_outside_2_3 = [v for v in solution if (v.is_real and Or(v<2,v>3))]
>>> print(solution_outside_2_3)
[1, 4]
对于 solveset()，在函数调用中通过设置一个域来限制输出的定义域。
>>> from sympy import S, solveset
>>> from sympy.abc import x
>>> solveset(x**4 - 256, x, domain=S.Reals)
{-4, 4}
或者通过将返回的解限制为任意集合，包括一个区间：
>>> from sympy import Interval, pi, sin, solveset
>>> from sympy.abc import x
>>> solveset(sin(x), x, Interval(-pi, pi))
{0, -pi, pi}
如果你将解限制在没有解的域中，solveset() 将返回空集合，EmptySet。
>>> from sympy import solveset, S
>>> from sympy.abc import x
>>> solveset(x**2 + 1, x, domain=S.Reals)
EmptySet
显式地表示可能解的无限集合。
solveset() 可以表示可能解的无限集合，并以标准数学符号表示，例如对于每个整数值的 (n)，满足 (\sin(x) = 0) 的 (x = n * \pi)：
>>> from sympy import pprint, sin, solveset
>>> from sympy.abc import x
>>> solution = solveset(sin(x), x)
>>> pprint(solution)
{2*n*pi | n in Integers} U {2*n*pi + pi | n in Integers}
然而，solve() 只会返回有限数量的解：
>>> from sympy import sin, solve
>>> from sympy.calculus.util import periodicity
>>> from sympy.abc import x
>>> f = sin(x)
>>> solve(f, x)
[0, pi]
>>> periodicity(f, x)
2*pi
solve() 尝试返回足够多的解，以便通过添加方程的周期性（此处为( 2\pi )的整数倍）生成所有（无穷多个）解。
使用解结果
将solve()的解代入表达式中
您可以将solve()的解代入表达式中。
一个常见的用例是找到函数( f )的临界点和值。在临界点，Derivative等于零（或未定义）。然后，您可以通过将临界点代入函数中使用subs()来获取这些临界点的函数值。您还可以通过将值代入二阶导数表达式来判断临界点是否为最大值或最小值：负值表示最大值，正值表示最小值。
>>> from sympy.abc import x
>>> from sympy import solve, diff
>>> f = x**3 + x**2 - x
>>> derivative = diff(f, x)
>>> critical_points = solve(derivative, x, dict=True)
>>> print(critical_points)
[{x: -1}, {x: 1/3}]
>>> point1, point2 = critical_points
>>> print(f.subs(point1))
1
>>> print(f.subs(point2))
-5/27
>>> curvature = diff(f, x, 2)
>>> print(curvature.subs(point1))
-4
>>> print(curvature.subs(point2))
4
solveset() 解集可能无法通过编程方式查询。
如果solveset()返回一个有限集（类FiniteSet)，您可以遍历解：
>>> from sympy import solveset
>>> from sympy.abc import x, y
>>> solution_set = solveset(x**2 - y, x)
>>> print(solution_set)
{-sqrt(y), sqrt(y)}
>>> solution_list = list(solution_set)
>>> print(solution_list)
[sqrt(y), -sqrt(y)]
然而，对于更复杂的结果，可能无法列出所有解。
>>> from sympy import S, solveset, symbols
>>> x, y = symbols('x, y')
>>> solution_set = solveset(x**2 - y, x, domain=S.Reals)
>>> print(solution_set)
Intersection({-sqrt(y), sqrt(y)}, Reals)
>>> list(solution_set)
Traceback (most recent call last):
...
TypeError: The computation had not completed because of the undecidable set
membership is found in every candidates.
在这种情况下，这是因为如果( y )为负数，其平方根将是虚数而不是实数，因此超出了解集的声明域。通过声明( y )为实数且为正，SymPy 可以确定其平方根为实数，从而解决解集与实数集之间的交集：
>>> from sympy import S, Symbol, solveset
>>> x = Symbol('x')
>>> y = Symbol('y', real=True, positive=True)
>>> solution_set = solveset(x**2 - y, x, domain=S.Reals)
>>> print(solution_set)
{-sqrt(y), sqrt(y)}
>>> list(solution_set)
[sqrt(y), -sqrt(y)]
或者，您可以从解集中提取集合，使用args，然后从包含符号解的集合中创建列表：
>>> from sympy import S, solveset, symbols
>>> x, y = symbols('x, y')
>>> solution_set = solveset(x**2 - y, x, domain=S.Reals)
>>> print(solution_set)
Intersection({-sqrt(y), sqrt(y)}, Reals)
>>> solution_set_args = solution_set.args
>>> print(solution_set.args)
(Reals, {-sqrt(y), sqrt(y)})
>>> list(solution_set_args[1])
[sqrt(y), -sqrt(y)]
可以加快solve()的选项
参考 solving guidance。
并非所有方程都能求解
没有封闭形式解的方程
有些方程没有闭式解，此时 SymPy 可能返回一个空集或出现错误。例如，下面的超越方程没有闭式解：
>>> from sympy import cos, solve
>>> from sympy.abc import x
>>> solve(cos(x) - x, x, dict=True)
Traceback (most recent call last):
...
NotImplementedError: multiple generators [x, cos(x)]
No algorithms are implemented to solve equation -x + cos(x)
有闭式解的方程，而 SymPy 无法解决
可能也存在一个代数解决方案来解决你的方程，但 SymPy 尚未实现适当的算法。如果发生这种情况，或者当 SymPy 返回一个空集或列表时（表示 SymPy 中存在 bug），请在邮件列表上发布，或在SymPy 的 GitHub 页面上开一个 issue。在问题解决之前，你可以数值解你的方程。
报告 Bug
如果你发现解决函数存在 bug，请在SymPy 邮件列表上发布问题。在问题解决之前，你可以使用考虑的备选方案中列出的其他方法。
代数方法求解方程组
原文：
docs.sympy.org/latest/guides/solving/solve-system-of-equations-algebraically.html
使用 SymPy 代数方法求解线性或非线性方程组。例如，对于解 (x² + y = 2z, y = -4z) 求解 x 和 y（假设 z 是常数或参数）得到 ({(x = -\sqrt{6z}, y = -4z),) ({(x = \sqrt{6z}, y = -4z)}})。
考虑的替代方案
- 一些方程组无法通过代数方法（无论是完全还是通过 SymPy）求解，因此您可能需要通过数值方法 数值求解您的方程组，而不是使用
nsolve()。
代数方法解决方程组的示例
无论您的方程是线性还是非线性，您都可以使用 solve():
代数方法求解线性方程组
>>> from sympy import solve
>>> from sympy.abc import x, y, z
>>> solve([x + y - 2*z, y + 4*z], [x, y], dict=True)
[{x: 6*z, y: -4*z}]
代数方法求解非线性方程组
>>> from sympy import solve
>>> from sympy.abc import x, y, z
>>> solve([x**2 + y - 2*z, y + 4*z], x, y, dict=True)
[{x: -sqrt(6)*sqrt(z), y: -4*z}, {x: sqrt(6)*sqrt(z), y: -4*z}]
指南
参考 在函数调用中包含要解决的变量 和 确保 solve() 的一致格式。
下面有两种方法来包含解决方案结果：字典 或 集合。字典更容易通过编程方式进行查询，因此如果需要使用代码提取解决方案，我们建议使用字典方法。
求解并使用结果作为一个字典
给出作为字典的解决方案
您可以为一些变量（例如，(x) 和 (y)) 解决一组方程，将另一个符号作为常数或参数（例如，(z))。您可以将要解决的变量指定为多个单独的参数，或作为一个列表（或元组）：
>>> from sympy import solve
>>> from sympy.abc import x, y, z
>>> equations = [x**2 + y - 2*z, y + 4*z]
>>> solutions = solve(equations, x, y, dict=True)
>>> solutions
[{x: -sqrt(6)*sqrt(z), y: -4*z}, {x: sqrt(6)*sqrt(z), y: -4*z}]
使用字典给出的解决方案
然后，您可以通过索引（用方括号指定）解的编号，然后是符号来提取解决方案。例如 solutions[0][x] 给出第一个解的 x 的结果：
>>> solutions[0][x]
-sqrt(6)*sqrt(z)
>>> solutions[0][y]
-4*z
求解结果为一个集合
若要获取符号列表和解集，请使用 set=True 而不是 dict=True：
from sympy import solve
from sympy.abc import x, y, z
solve([x**2 + y - 2*z, y + 4*z], [x, y], set=True)
([x, y], {(-sqrt(6)*sqrt(z), -4*z), (sqrt(6)*sqrt(z), -4*z)})
加快 solve() 的选项
参考 加快 solve() 的选项。
并非所有的方程组都可以求解
无解的方程组
一些方程组没有解。例如，以下两个方程组没有解，因为它们归结为 1 == 0，所以 SymPy 返回一个空列表：
>>> from sympy import solve
>>> from sympy.abc import x, y
>>> solve([x + y - 1, x + y], [x, y], dict=True)
[]
from sympy import solve
from sympy.abc import x, y, z
solve([x + y - (z + 1), x + y - z)], [x, y], dict=True)
[]
下面的系统简化为 (z = 2z)，因此没有通解，但如果 (z=0)，则可能满足。请注意，solve() 不会假定 (z=0)，即使这是使方程组一致的唯一值，因为 (z) 是一个参数而不是未知数。也就是说，solve() 不会像处理未知数一样处理 (z)，因为它不在指定为未知数解的符号列表中（[x, y]），所有这些符号都像具有任意值的参数一样处理。一个符号是变量还是参数的区分只有在使用symbols()（或从abc导入）创建符号时才能确定。在创建符号时，没有这样的区别。
>>> from sympy import solve
>>> from sympy.abc import x, y, z
>>> solve([x + y - z, x + y - 2*z], [x, y], dict=True)
[]
下面的系统是过约束的，意味着方程（三个）比要解的未知数（两个，即 (x) 和 (y)）更多。它没有解：
>>> from sympy import solve
>>> from sympy.abc import x, y, z
>>> solve([x + y - z, x - (z + 1), 2*x - y], [x, y], dict=True)
[]
注意，一些过约束系统确实有解（例如，如果一个方程是其他方程的线性组合），在这种情况下 SymPy 可以解决这种过约束系统。
没有封闭形式解的方程组
某些方程组无法通过代数方法求解，例如包含超越方程的方程组：
>>> from sympy import cos, solve
>>> from sympy.abc import x, y, z
>>> solve([x - y, cos(x) - y], [x, y], dict=True)
Traceback (most recent call last):
...
NotImplementedError: could not solve -y + cos(y)
所以你可以使用nsolve()来找到数值解：
>>> from sympy import cos, nsolve
>>> from sympy.abc import x, y, z
>>> nsolve([x - y, cos(x) - y], [x, y], [1,1])
Matrix([
[0.739085133215161],
[0.739085133215161]])
有封闭形式解但 SymPy 无法解决的方程：
也可能是你的方程有代数解，但 SymPy 尚未实现适当的算法。如果 SymPy 返回一个空集或列表，而你知道存在封闭形式解（表明 SymPy 存在错误），请在邮件列表上发布问题，或在SymPy 的 GitHub 页面上开一个问题。在问题解决之前，可以使用考虑的替代方法中列出的其他方法。
报告 Bug
如果你发现solve()存在 bug，请在SymPy 邮件列表上发布问题。在问题解决之前，可以使用考虑的替代方法中列出的其他方法。
数值求解一个或一组方程
原文：
docs.sympy.org/latest/guides/solving/solve-numerically.html
使用 SymPy 来数值求解一个或多个方程组。例如，数值解 (\cos(x) = x ) 返回 ( x \approx 0.739085133215161)。
如果需要：
-
您只需要一个数值解，而不是符号解
-
没有可用的闭式解或者解法过于复杂；参考何时可能更喜欢数值解
solve() 和 solveset() 不会尝试找到数值解，只会找到数学上精确的符号解。因此，如果您需要数值解，请使用 nsolve()。
SymPy 是为符号数学设计的。如果您不需要进行符号操作，则对于数值运算，可以使用另一个免费开源的软件包，如 NumPy 或 SciPy，它们速度更快，适用于数组，并且实现了更多的算法。使用 SymPy（或其依赖项 mpmath）进行数值计算的主要原因是：
-
在 SymPy 中进行符号计算的上下文中进行简单的数值计算
-
如果您需要任意精度功能以获得比 float64 更多位数的精度。
考虑的替代方案
-
SciPy 的
scipy.optimize.fsolve()可以解决一组（非线性）方程 -
NumPy 的
numpy.linalg.solve()可以解决一组线性标量方程 -
mpmath 的
findroot()，它被nsolve()调用并可以传递参数给它
数值求解方程的示例
这是一个数值解决一个方程的示例：
>>> from sympy import cos, nsolve, Symbol
>>> x = Symbol('x')
>>> nsolve(cos(x) - x, x, 1)
0.739085133215161
指导
支持过度确定的方程组。
找到一个实函数的复根
要解决实函数的复根，请指定一个非实数（纯虚数或复数）的初始点：
>>> from sympy import nsolve
>>> from sympy.abc import x
>>> nsolve(x**2 + 2, 1) # Real initial point returns no root
Traceback (most recent call last):
...
ValueError: Could not find root within given tolerance. (4.18466446988997098217 > 2.16840434497100886801e-19)
Try another starting point or tweak arguments.
>>> from sympy import I
>>> nsolve(x**2 + 2, I) # Imaginary initial point returns a complex root
1.4142135623731*I
>>> nsolve(x**2 + 2, 1 + I) # Complex initial point returns a complex root
1.4142135623731*I
确保找到的根在给定区间内
不保证nsolve()会找到距离初始点最近的根。在这里，即使根-1更接近初始点-0.1，nsolve()也找到了根1：
>>> from sympy import nsolve
>>> from sympy.abc import x
>>> nsolve(x**2 - 1, -0.1)
1.00000000000000
您可以通过指定一个元组中的区间，并使用 solver='bisect' 来确保找到的根位于给定区间内（如果存在这样的根）。在这里，指定区间 (-10, 0) 确保找到根-1：
>>> from sympy import nsolve
>>> from sympy.abc import x
>>> nsolve(x**2 - 1, (-10, 0), solver='bisect')
-1.00000000000000
数值解决方程组
要解决多维函数系统，请提供一个元组
-
函数
(f1, f2) -
变量解为
(x1, x2) -
起始值
(-1, 1)
>>> from sympy import Symbol, nsolve
>>> x1 = Symbol('x1')
>>> x2 = Symbol('x2')
>>> f1 = 3 * x1**2 - 2 * x2**2 - 1
>>> f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
>>> print(nsolve((f1, f2), (x1, x2), (-1, 1)))
Matrix([[-1.19287309935246], [1.27844411169911]])
提高解的精度
您可以使用 prec 来增加解的精度：
>>> from sympy import Symbol, nsolve
>>> x1 = Symbol('x1')
>>> x2 = Symbol('x2')
>>> f1 = 3 * x1**2 - 2 * x2**2 - 1
>>> f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
>>> print(nsolve((f1, f2), (x1, x2), (-1, 1), prec=25))
Matrix([[-1.192873099352460791205211], [1.278444111699106966687122]])
创建可以使用 SciPy 求解的函数
如上所述，SymPy 专注于符号计算，不适用于数值计算。如果需要频繁调用数值求解器，则最好使用专为数值计算优化的求解器，如 SciPy 的 root_scalar()。推荐的工作流程是：
-
使用 SymPy 生成（通过符号化简或解决方程）数学表达式
-
使用
lambdify()将其转换为 lambda 函数 -
使用类似 SciPy 的数值库生成数值解
>>> from sympy import simplify, cos, sin, lambdify
>>> from sympy.abc import x, y
>>> from scipy.optimize import root_scalar
>>> expr = cos(x * (x + x**2)/(x*sin(y)**2 + x*cos(y)**2 + x))
>>> simplify(expr) # 1\. symbolically simplify expression
cos(x*(x + 1)/2)
>>> lam_f = lambdify(x, cos(x*(x + 1)/2)) # 2\. lambdify
>>> sol = root_scalar(lam_f, bracket=[0, 2]) # 3\. numerically solve using SciPy
>>> sol.root
1.3416277185114782
使用解的结果
将结果代入表达式中
最佳做法是使用 evalf() 将数值值替换为表达式。以下代码示例表明，数值值并非精确的根，因为将其代回表达式会产生一个与零略有不同的结果：
>>> from sympy import cos, nsolve, Symbol
>>> x = Symbol('x')
>>> f = cos(x) - x
>>> x_value = nsolve(f, x, 1); x_value
0.739085133215161
>>> f.evalf(subs={x: x_value})
-5.12757857962640e-17
使用 subs 可能会由于精度误差而得到错误的结果，在这里将 -5.12757857962640e-17 有效地舍入为零：
>>> f.subs(x, x_value)
0
在替换值时，可以将一些符号留作变量：
>>> from sympy import cos, nsolve, Symbol
>>> x = Symbol('x')
>>> f = cos(x) - x
>>> x_value = nsolve(f, x, 1); x_value
0.739085133215161
>>> y = Symbol('y')
>>> z = Symbol('z')
>>> g = x * y**2
>>> values = {x: x_value, y: 1}
>>> (x + y - z).evalf(subs=values)
1.73908513321516 - z
并非所有方程都可以解决
nsolve() 是一个数值求解函数，因此它经常可以为无法代数求解的方程提供解决方案。
没有解的方程
一些方程无解，这种情况下 SymPy 可能会返回错误。例如，方程 (e^x = 0)（在 SymPy 中为 exp(x)）无解：
>>> from sympy import nsolve, exp
>>> from sympy.abc import x
>>> nsolve(exp(x), x, 1, prec=20)
Traceback (most recent call last):
...
ValueError: Could not find root within given tolerance. (5.4877893607115270300540019e-18 > 1.6543612251060553497428174e-24)
Try another starting point or tweak arguments.
报告错误
如果你在使用nsolve()时发现了 bug，请在SymPy 邮件列表上发布问题。在问题解决之前，你可以考虑使用备选方法中列出的其他方法。
代数解一个常微分方程（ODE）
使用 SymPy 代数解一个常微分方程（ODE）。例如，解(y''(x) + 9y(x)=0 )得到( y(x)=C_{1} \sin(3x)+ C_{2} \cos(3x))。
可供考虑的替代方案
- 要数值解一个 ODE 系统，可以使用 SciPy ODE solver ，如
solve_ivp。你也可以使用 SymPy 创建然后使用 SciPy 的solve_ivp数值求解一个 ODE，具体可见下文的在 SciPy 中数值解 ODE。
解一个常微分方程（ODE）
这里是使用dsolve() 代数解上述 ODE 的示例。你可以使用checkodesol() 来验证解是否正确。
>>> from sympy import Function, dsolve, Derivative, checkodesol
>>> from sympy.abc import x
>>> y = Function('y')
>>> # Solve the ODE
>>> result = dsolve(Derivative(y(x), x, x) + 9*y(x), y(x))
>>> result
Eq(y(x), C1*sin(3*x) + C2*cos(3*x))
>>> # Check that the solution is correct
>>> checkodesol(Derivative(y(x), x, x) + 9*y(x), result)
(True, 0)
checkodesol() 的输出是一个元组，其中第一项是布尔值，指示将解代入 ODE 是否结果为0，表示解是正确的。
指导
定义导数
表达函数的导数有许多方式。对于未定义的函数，Derivative 和 diff() 都表示未定义的导数。因此，所有以下的 ypp (“y prime prime”) 都代表(y'')，即函数(y(x))关于(x)的二阶导数：
ypp = y(x).diff(x, x)
ypp = y(x).diff(x, 2)
ypp = y(x).diff((x, 2))
ypp = diff(y(x), x, x)
ypp = diff(y(x), x, 2)
ypp = Derivative(y(x), x, x)
ypp = Derivative(y(x), x, 2)
ypp = Derivative(Derivative(y(x), x), x)
ypp = diff(diff(y(x), x), x)
yp = y(x).diff(x)
ypp = yp.diff(x)
我们建议将待解函数作为dsolve() 的第二个参数进行指定。请注意，这必须是一个函数而不是一个变量（符号）。如果你指定了一个变量（(x)）而不是一个函数（(f(x))），SymPy 将会报错：
>>> dsolve(Derivative(y(x), x, x) + 9*y(x), x)
Traceback (most recent call last):
...
ValueError: dsolve() and classify_ode() only work with functions of one variable, not x
同样，你必须指定函数的参数：(y(x))，而不仅仅是(y)。
定义 ODE 的选项
你可以通过两种方式定义待解的函数。对于你选择的初始条件指定语法取决于你的选择。
选项 1：定义一个不包括其自变量的函数
你可以定义一个不包括其自变量的函数：
>>> from sympy import symbols, Eq, Function, dsolve
>>> f, g = symbols("f g", cls=Function)
>>> x = symbols("x")
>>> eqs = [Eq(f(x).diff(x), g(x)), Eq(g(x).diff(x), f(x))]
>>> dsolve(eqs, [f(x), g(x)])
[Eq(f(x), -C1*exp(-x) + C2*exp(x)), Eq(g(x), C1*exp(-x) + C2*exp(x))]
请注意，作为dsolve() 的第二个参数，你需要提供待解的函数列表，如此处 [f(x), g(x)]。
指定初始条件或边界条件
如果您的微分方程具有初始条件或边界条件，请使用 dsolve() 的可选参数 ics 指定它们。初始条件和边界条件以相同方式处理（尽管参数称为 ics）。应以 {f(x0): y0, f(x).diff(x).subs(x, x1): y1} 形式给出，例如在 (x = x_{0}) 处 (f(x)) 的值是 (y_{0})。对于幂级数解，如果未指定初始条件，则假定 (f(0)) 为 (C_{0})，并且关于 (0) 计算幂级数解。
这里是一个设置函数初始值的例子，即 (f(0) = 1) 和 (g(2) = 3)：
>>> from sympy import symbols, Eq, Function, dsolve
>>> f, g = symbols("f g", cls=Function)
>>> x = symbols("x")
>>> eqs = [Eq(f(x).diff(x), g(x)), Eq(g(x).diff(x), f(x))]
>>> dsolve(eqs, [f(x), g(x)])
[Eq(f(x), -C1*exp(-x) + C2*exp(x)), Eq(g(x), C1*exp(-x) + C2*exp(x))]
>>> dsolve(eqs, [f(x), g(x)], ics={f(0): 1, g(2): 3})
[Eq(f(x), (1 + 3*exp(2))*exp(x)/(1 + exp(4)) - (-exp(4) + 3*exp(2))*exp(-x)/(1 + exp(4))), Eq(g(x), (1 + 3*exp(2))*exp(x)/(1 + exp(4)) + (-exp(4) + 3*exp(2))*exp(-x)/(1 + exp(4)))]
这里是设置函数导数初始值的例子，即 (f'(1) = 2)：
>>> eqn = Eq(f(x).diff(x), f(x))
>>> dsolve(eqn, f(x), ics={f(x).diff(x).subs(x, 1): 2})
Eq(f(x), 2*exp(-1)*exp(x))
选项 2：定义一个独立变量的函数
您可能更喜欢指定一个函数（例如 (y)) 的独立变量（例如 (t))，这样 y 就表示 y(t)：
>>> from sympy import symbols, Function, dsolve
>>> t = symbols('t')
>>> y = Function('y')(t)
>>> y
y(t)
>>> yp = y.diff(t)
>>> ypp = yp.diff(t)
>>> eq = ypp + 2*yp + y
>>> eq
y(t) + 2*Derivative(y(t), t) + Derivative(y(t), (t, 2))
>>> dsolve(eq, y)
Eq(y(t), (C1 + C2*t)*exp(-t))
使用此约定，dsolve() 的第二个参数 y 表示 y(t)，因此 SymPy 将其识别为要解的有效函数。
指定初始条件或边界条件
使用该语法，您可以通过使用 subs() 将独立变量的值替换到函数 (y) 中，因为函数 (y) 已经将其独立变量作为参数 (t)：
>>> dsolve(eq, y, ics={y.subs(t, 0): 0})
Eq(y(t), C2*t*exp(-t))
注意复制和粘贴结果
如果您选择定义一个独立变量的函数，请注意复制结果并粘贴到后续代码中可能会导致错误，因为 x 已经定义为 y(t)，所以如果您粘贴 y(t)，它会被解释为 y(t)(t)：
>>> dsolve(y(t).diff(y), y)
Traceback (most recent call last):
...
TypeError: 'y' object is not callable
因此，请记住不要包含独立变量调用 (t)：
>>> dsolve(y.diff(t), y)
Eq(y(t), C1)
使用解决方案结果
不同于其他求解函数，dsolve() 返回一个以如下格式的 Equality（方程）：Eq(y(x), C1*sin(3*x) + C2*cos(3*x))，这等同于数学符号 (y(x) = C_1 \sin(3x) + C_2 \cos(3x))。
提取单个解和函数的结果
您可以从 Equality 中使用右侧属性 rhs 提取结果：
>>> from sympy import Function, dsolve, Derivative
>>> from sympy.abc import x
>>> y = Function('y')
>>> result = dsolve(Derivative(y(x), x, x) + 9*y(x), y(x))
>>> result
Eq(y(x), C1*sin(3*x) + C2*cos(3*x))
>>> result.rhs
C1*sin(3*x) + C2*cos(3*x)
有些常微分方程不能显式求解，只能隐式求解
上述常微分方程可以显式求解，特别是 (y(x)) 可以用 (x) 的函数表示。然而，有些常微分方程不能显式求解，例如：
>>> from sympy import dsolve, exp, symbols, Function
>>> f = symbols("f", cls=Function)
>>> x = symbols("x")
>>> dsolve(f(x).diff(x) + exp(-f(x))*f(x))
Eq(Ei(f(x)), C1 - x)
这不直接给出了 (f(x)) 的表达式。相反，dsolve() 将一个解表达为 (g(f(x)))，其中 (g) 是Ei，经典指数积分函数。Ei 没有已知的闭合形式逆运算，所以一个解不能明确地表达为 (f(x)) 等于 (x) 的函数。相反，dsolve 返回一个隐式解。
当dsolve返回一个隐式解时，提取返回的等式的右侧将不会给出一个明确的表达式，用于要解的函数，这里是(f(x))。因此，在提取要解的函数的表达式之前，检查dsolve能否明确为该函数求解。
提取多个函数-解对的结果
如果您正在解决一个具有多个未知函数的方程组，dsolve()的输出形式取决于是否有一个或多个解。
如果存在一个解集
如果一个多元函数方程组只有一个解集，dsolve()将返回一个非嵌套的包含一个等式的列表。您可以使用单个循环或推导式提取解表达式：
>>> from sympy import symbols, Eq, Function, dsolve
>>> y, z = symbols("y z", cls=Function)
>>> x = symbols("x")
>>> eqs_one_soln_set = [Eq(y(x).diff(x), z(x)**2), Eq(z(x).diff(x), z(x))]
>>> solutions_one_soln_set = dsolve(eqs_one_soln_set, [y(x), z(x)])
>>> solutions_one_soln_set
[Eq(y(x), C1 + C2**2*exp(2*x)/2), Eq(z(x), C2*exp(x))]
>>> # Loop through list approach
>>> solution_one_soln_set_dict = {}
>>> for fn in solutions_one_soln_set:
... solution_one_soln_set_dict.update({fn.lhs: fn.rhs})
>>> solution_one_soln_set_dict
{y(x): C1 + C2**2*exp(2*x)/2, z(x): C2*exp(x)}
>>> # List comprehension approach
>>> solution_one_soln_set_dict = {fn.lhs:fn.rhs for fn in solutions_one_soln_set}
>>> solution_one_soln_set_dict
{y(x): C1 + C2**2*exp(2*x)/2, z(x): C2*exp(x)}
>>> # Extract expression for y(x)
>>> solution_one_soln_set_dict[y(x)]
C1 + C2**2*exp(2*x)/2
如果存在多个解集
如果一个多元函数方程组有多个解集，dsolve() 将返回一个嵌套的等式列表，外部列表表示每个解，内部列表表示每个函数。虽然您可以通过指定每个函数的索引来提取结果，但我们建议一种对函数排序具有鲁棒性的方法。以下将每个解转换为字典，以便您可以轻松提取所需函数的结果。它使用标准的 Python 技术，如循环或推导式，以嵌套的方式。
>>> from sympy import symbols, Eq, Function, dsolve
>>> y, z = symbols("y z", cls=Function)
>>> x = symbols("x")
>>> eqs = [Eq(y(x).diff(x)**2, z(x)**2), Eq(z(x).diff(x), z(x))]
>>> solutions = dsolve(eqs, [y(x), z(x)])
>>> solutions
[[Eq(y(x), C1 - C2*exp(x)), Eq(z(x), C2*exp(x))], [Eq(y(x), C1 + C2*exp(x)), Eq(z(x), C2*exp(x))]]
>>> # Nested list approach
>>> solutions_list = []
>>> for solution in solutions:
... solution_dict = {}
... for fn in solution:
... solution_dict.update({fn.lhs: fn.rhs})
... solutions_list.append(solution_dict)
>>> solutions_list
[{y(x): C1 - C2*exp(x), z(x): C2*exp(x)}, {y(x): C1 + C2*exp(x), z(x): C2*exp(x)}]
>>> # Nested comprehension approach
>>> solutions_list = [{fn.lhs:fn.rhs for fn in solution} for solution in solutions]
>>> solutions_list
[{y(x): C1 - C2*exp(x), z(x): C2*exp(x)}, {y(x): C1 + C2*exp(x), z(x): C2*exp(x)}]
>>> # Extract expression for y(x)
>>> solutions_list[0][y(x)]
C1 - C2*exp(x)
处理任意常数
您可以操纵由dsolve()自动生成的C1、C2和C3等任意常数，方法是将它们创建为符号。例如，如果您想为任意常数分配值，可以将它们创建为符号，然后使用subs()替换它们的值：
>>> from sympy import Function, dsolve, Derivative, symbols, pi
>>> y = Function('y')
>>> x, C1, C2 = symbols("x, C1, C2")
>>> result = dsolve(Derivative(y(x), x, x) + 9*y(x), y(x)).rhs
>>> result
C1*sin(3*x) + C2*cos(3*x)
>>> result.subs({C1: 7, C2: pi})
7*sin(3*x) + pi*cos(3*x)
在 SciPy 中数值求解 ODE
利用SciPy快速数值 ODE 求解的一种常见工作流程是
-
在 SymPy 中设置一个 ODE
-
使用
lambdify()将其转换为数值函数。 -
通过使用 SciPy 的
solve_ivp数值积分 ODE 来解决初值问题来解决初始值问题。
这里是关于化学动力学领域的示例，其中非线性常微分方程采用以下形式：
[\begin{split} r_f = & k_f y_0(t)² y_1(t) \ r_b = & k_b y_2(t)² \ \frac{d y_0(t)}{dt} = & 2(r_b - r_f) \ \frac{d y_1(t)}{dt} = & r_b - r_f \ \frac{d y_2(t)}{dt} = & 2(r_f - r_b) \end{split}]
和
[\begin{split}\vec{y}(t) = \begin{bmatrix} y_0(t) \ y_1(t) \ y_2(t) \end{bmatrix} \end{split}]
>>> from sympy import symbols, lambdify
>>> import numpy as np
>>> import scipy.integrate
>>> import matplotlib.pyplot as plt
>>> # Create symbols y0, y1, and y2
>>> y = symbols('y:3')
>>> kf, kb = symbols('kf kb')
>>> rf = kf * y[0]**2 * y[1]
>>> rb = kb * y[2]**2
>>> # Derivative of the function y(t); values for the three chemical species
>>> # for input values y, kf, and kb
>>> ydot = [2*(rb - rf), rb - rf, 2*(rf - rb)]
>>> ydot
[2*kb*y2**2 - 2*kf*y0**2*y1, kb*y2**2 - kf*y0**2*y1, -2*kb*y2**2 + 2*kf*y0**2*y1]
>>> t = symbols('t') # not used in this case
>>> # Convert the SymPy symbolic expression for ydot into a form that
>>> # SciPy can evaluate numerically, f
>>> f = lambdify((t, y, kf, kb), ydot)
>>> k_vals = np.array([0.42, 0.17]) # arbitrary in this case
>>> y0 = [1, 1, 0] # initial condition (initial values)
>>> t_eval = np.linspace(0, 10, 50) # evaluate integral from t = 0-10 for 50 points
>>> # Call SciPy's ODE initial value problem solver solve_ivp by passing it
>>> # the function f,
>>> # the interval of integration,
>>> # the initial state, and
>>> # the arguments to pass to the function f
>>> solution = scipy.integrate.solve_ivp(f, (0, 10), y0, t_eval=t_eval, args=k_vals)
>>> # Extract the y (concentration) values from SciPy solution result
>>> y = solution.y
>>> # Plot the result graphically using matplotlib
>>> plt.plot(t_eval, y.T)
>>> # Add title, legend, and axis labels to the plot
>>> plt.title('Chemical Kinetics')
>>> plt.legend(['NO', 'Br$_2$', 'NOBr'], shadow=True)
>>> plt.xlabel('time')
>>> plt.ylabel('concentration')
>>> # Finally, display the annotated plot
>>> plt.show()
(png, hires.png, pdf)
SciPy 的solve_ivp返回一个结果，其中包含每个化学物种对应于时间点t_eval的y（数值函数结果，这里是浓度）值。
普通微分方程求解提示
返回未评估的积分
默认情况下，dsolve()尝试评估它生成的积分以解决您的普通微分方程。您可以通过使用以_Integral结尾的提示函数来禁用积分的评估，例如separable_Integral。这是有用的，因为integrate()是一个昂贵的例程。由于难以或无法积分，SymPy 可能会挂起（似乎永远无法完成操作），因此使用_Integral提示至少会返回一个（未积分的）结果，您可以随后考虑。禁用积分的最简单方法是使用all_Integral提示，因为您不需要知道要提供哪种提示：对于具有相应的_Integral提示的任何提示，all_Integral只返回_Integral提示。
选择特定的求解器
您可能希望选择特定的求解器，有几个原因：
-
教育目的：例如，如果您正在学习某种特定的解 ODE 方法，并希望获得完全匹配该方法的结果
-
结果形式：有时候一个常微分方程可以由许多不同的求解器求解，它们可以返回不同的结果。尽管它们在数学上是等价的，但任意常数可能不同。
dsolve()默认情况下会尝试首先使用“最佳”求解器，这些求解器最有可能产生最有用的输出，但这不是一个完美的启发式。例如，“最佳”求解器可能生成一个包含 SymPy 无法解决的积分的结果，但另一个求解器可能生成一个 SymPy 可以解决的不同积分。因此，如果解决方案不符合您的要求，您可以尝试其他提示，以查看它们是否提供更好的结果。
并非所有方程都可以解决
没有解的方程
并非所有的微分方程都可以解决，例如：
>>> from sympy import Function, dsolve, Derivative, symbols
>>> y = Function('y')
>>> x, C1, C2 = symbols("x, C1, C2")
>>> dsolve(Derivative(y(x), x, 3) - (y(x)**2), y(x)).rhs
Traceback (most recent call last):
...
NotImplementedError: solve: Cannot solve -y(x)**2 + Derivative(y(x), (x, 3))
没有封闭形式解的方程
如上所述，有些常微分方程只能隐式求解。
有些微分方程组没有封闭形式的解决方案，因为它们是混沌的，例如Lorenz system或由以下这两个微分方程描述的双摆（从ScienceWorld简化而来）：
[ 2 \theta_1''(t) + \theta_2''(t) \cos(\theta_1-\theta_2) + \theta_2'²(t) \sin(\theta_1 - \theta_2) + 2g \sin(\theta_1) = 0 ][ \theta_2''(t) + \theta_1''(t) \cos(\theta_1-\theta_2) - \theta_1'²(t) \sin(\theta_1 - \theta_2) + g \sin(\theta_2) = 0 ]
>>> from sympy import symbols, Function, cos, sin, dsolve
>>> theta1, theta2 = symbols('theta1 theta2', cls=Function)
>>> g, t = symbols('g t')
>>> eq1 = 2*theta1(t).diff(t, t) + theta2(t).diff(t, t)*cos(theta1(t) - theta2(t)) + theta2(t).diff(t)**2*sin(theta1(t) - theta2(t)) + 2*g*sin(theta1(t))
>>> eq2 = theta2(t).diff(t, t) + theta1(t).diff(t, t)*cos(theta1(t) - theta2(t)) - theta1(t).diff(t)**2*sin(theta1(t) - theta2(t)) + g*sin(theta2(t))
>>> dsolve([eq1, eq2], [theta1(t), theta2(t)])
Traceback (most recent call last):
...
NotImplementedError
对于这种情况，您可以如考虑的替代方案中提到的那样通过数值方法来解方程。
报告错误
如果您在dsolve()中发现了一个 bug，请在SymPy 邮件列表上发布问题。在问题解决之前，您可以使用考虑的替代方法中列出的其他方法。
代数或数值上找到多项式的根
原文：
docs.sympy.org/latest/guides/solving/find-roots-polynomial.html
使用 SymPy 代数上找到一元多项式的根。例如，对于 (ax² + bx + c) 找到 (x) 的根为 (x = \frac{-b\pm\sqrt{b² - 4ac}}{2a})。
考虑的替代方案
代数上找到多项式根的例子
这里是一个代数上找到多项式根的例子：
>>> from sympy import roots
>>> from sympy.abc import x, a, b, c
>>> roots(a*x**2 + b*x + c, x)
{-b/(2*a) - sqrt(-4*a*c + b**2)/(2*a): 1,
-b/(2*a) + sqrt(-4*a*c + b**2)/(2*a): 1}
此示例重现了二次公式。
找到多项式根的函数
有几个函数可以用来找到多项式的根：
-
solve()是一个通用的解函数，可以找到根，但效率低于all_roots()，并且是此列表中唯一不传达根的重数的函数；solve()也适用于非多项式方程和非多项式方程组 -
roots()计算一元多项式的符号根；对于大多数高次多项式（五次或更高次）将失败。 -
nroots()计算可以数值评估系数的任何多项式的数值近似根，无论系数是有理数还是无理数。 -
RootOf()可以精确表示任意大次数多项式的所有根，只要系数是有理数。RootOf()可以避免病态条件和返回虚假的复数部分，因为它使用基于隔离区间的更精确但更慢的数值算法。以下两个函数使用RootOf()，因此它们具有相同的属性：-
real_roots()可以精确找到任意大次数多项式的所有实根；因为它只找到实根，所以它可能比找到所有根的函数更有效。 -
all_roots()可以精确找到任意大次数多项式的所有根。
-
-
factor()将多项式因式分解为不可约多项式，并且可以显示根位于系数环中
每一个都将在本页上使用。
指导
参考 在函数调用中包含要解决的变量 和 使用精确值。
找到多项式的根
你可以通过多种方式代数地找到多项式的根。使用哪一种取决于你是否
-
想要代数或数值答案
-
想要每个根的重数（每个根是解的次数）。在表示为 ((x+2)²(x-3)) 的
expression下，根 -2 的重数为二，因为 (x+2) 被平方，而根 3 的重数为一，因为 (x-3) 没有指数。类似地，在symbolic表达式中，根 (-a) 的重数为二，根 (b) 的重数为一。
>>> from sympy import solve, roots, real_roots, factor, nroots, RootOf, expand
>>> from sympy import Poly
>>> expression = (x+2)**2 * (x-3)
>>> symbolic = (x+a)**2 * (x-b)
代数解决方案不考虑根重数。
你可以使用 SymPy 的标准 solve() 函数，尽管它不会返回根的重数：
>>> solve(expression, x, dict=True)
[{x: -2}, {x: 3}]
>>> solve(symbolic, x, dict=True)
[{x: -a}, {x: b}]
solve()首先尝试使用roots()；如果这不起作用，它将尝试使用all_roots()。对于三次（三次多项式）和四次（四次多项式），这意味着solve()将使用来自根的根式公式，而不是RootOf()即使 RootOf 是可能的。三次和四次公式通常给出在实际中无用的非常复杂的表达式。因此，您可能希望将solve()参数cubics或quartics设置为False以返回RootOf()结果：
>>> from sympy import solve
>>> from sympy.abc import x
>>> # By default, solve() uses the radical formula, yielding very complex terms
>>> solve(x**4 - x + 1, x)
[-sqrt(2/(3*(1/16 + sqrt(687)*I/144)**(1/3)) + 2*(1/16 + sqrt(687)*I/144)**(1/3))/2 - sqrt(-2*(1/16 + sqrt(687)*I/144)**(1/3) - 2/sqrt(2/(3*(1/16 + sqrt(687)*I/144)**(1/3)) + 2*(1/16 + sqrt(687)*I/144)**(1/3)) - 2/(3*(1/16 + sqrt(687)*I/144)**(1/3)))/2,
sqrt(2/(3*(1/16 + sqrt(687)*I/144)**(1/3)) + 2*(1/16 + sqrt(687)*I/144)**(1/3))/2 - sqrt(-2*(1/16 + sqrt(687)*I/144)**(1/3) + 2/sqrt(2/(3*(1/16 + sqrt(687)*I/144)**(1/3)) + 2*(1/16 + sqrt(687)*I/144)**(1/3)) - 2/(3*(1/16 + sqrt(687)*I/144)**(1/3)))/2,
sqrt(-2*(1/16 + sqrt(687)*I/144)**(1/3) - 2/sqrt(2/(3*(1/16 + sqrt(687)*I/144)**(1/3)) + 2*(1/16 + sqrt(687)*I/144)**(1/3)) - 2/(3*(1/16 + sqrt(687)*I/144)**(1/3)))/2 - sqrt(2/(3*(1/16 + sqrt(687)*I/144)**(1/3)) + 2*(1/16 + sqrt(687)*I/144)**(1/3))/2,
sqrt(-2*(1/16 + sqrt(687)*I/144)**(1/3) + 2/sqrt(2/(3*(1/16 + sqrt(687)*I/144)**(1/3)) + 2*(1/16 + sqrt(687)*I/144)**(1/3)) - 2/(3*(1/16 + sqrt(687)*I/144)**(1/3)))/2 + sqrt(2/(3*(1/16 + sqrt(687)*I/144)**(1/3)) + 2*(1/16 + sqrt(687)*I/144)**(1/3))/2]
>>> # If you set quartics=False, solve() uses RootOf()
>>> solve(x**4 - x + 1, x, quartics=False)
[CRootOf(x**4 - x + 1, 0),
CRootOf(x**4 - x + 1, 1),
CRootOf(x**4 - x + 1, 2),
CRootOf(x**4 - x + 1, 3)]
从标准数学符号中写出solve()的第一个根强调了它的复杂性：
[- \frac{\sqrt{\frac{2}{3 \sqrt[3]{\frac{1}{16} + \frac{\sqrt{687} i}{144}}} + 2 \sqrt[3]{\frac{1}{16} + \frac{\sqrt{687} i}{144}}}}{2} - \frac{\sqrt{- 2 \sqrt[3]{\frac{1}{16} + \frac{\sqrt{687} i}{144}} - \frac{2}{\sqrt{\frac{2}{3 \sqrt[3]{\frac{1}{16} + \frac{\sqrt{687} i}{144}}} + 2 \sqrt[3]{\frac{1}{16} + \frac{\sqrt{687} i}{144}}}} - \frac{2}{3 \sqrt[3]{\frac{1}{16} + \frac{\sqrt{687} i}{144}}}}}{2}]
此外，对于五次（五次多项式）或更高次多项式，没有一般的根式公式，因此它们的RootOf()表示可能是最佳选项。
参见代数解方程以了解更多关于使用solve()的信息。
代数解与根的重数
roots
roots()可以为具有符号系数的多项式的根给出显式表达式（即如果系数中有符号）。如果factor()没有揭示它们，则可能会失败。以下是roots()的示例：
>>> roots(expression, x)
{-2: 2, 3: 1}
>>> roots(symbolic, x)
{-a: 2, b: 1}
它以字典形式返回结果，其中键是根（例如，-2），值是该根的重数（例如，2）。
roots() 函数使用多种技术（因式分解、分解、根式公式）寻找根的表达式，如果可能的话返回根的根式表达式。当它能找到一些根的根式表达式时，它会返回它们及其重数。对于大多数高次多项式（五次或更高次），此函数将失败，因为它们没有根式解，并且不能保证它们根本上有闭合形式的解，这与阿贝尔-鲁菲尼定理的解释相符。
因式分解方程
另一种方法是使用factor()来因式分解多项式，它不直接给出根，但可以给出更简单的表达式：
>>> expression_expanded = expand(expression)
>>> expression_expanded
x**3 + x**2 - 8*x - 12
>>> factor(expression_expanded)
(x - 3)*(x + 2)**2
>>> symbolic_expanded = expand(symbolic)
>>> symbolic_expanded
-a**2*b + a**2*x - 2*a*b*x + 2*a*x**2 - b*x**2 + x**3
>>> factor(symbolic_expanded)
(a + x)**2*(-b + x)
factor() 也可以因式分解给定多项式环，这可以揭示根位于系数环中的信息。例如，如果多项式具有有理系数，则 factor() 将显示任何有理根。如果系数是涉及符号 (a) 的多项式，例如具有有理系数的多项式函数，则将显示与 (a) 有关的多项式函数的任何根。在此示例中，factor() 显示 (x = a²) 和 (x = -a³ - a) 是根：
>>> from sympy import expand, factor
>>> from sympy.abc import x, a
>>> p = expand((x - a**2)*(x + a + a**3))
>>> p
-a**5 + a**3*x - a**3 - a**2*x + a*x + x**2
>>> factor(p)
(-a**2 + x)*(a**3 + a + x)
精确数值解与根重数
real_roots
如果多项式的根是实数，使用real_roots() 确保只返回实数根（不包括复数或虚数）。
>>> from sympy import real_roots
>>> from sympy.abc import x
>>> cubed = x**3 - 1
>>> # roots() returns real and complex roots
>>> roots(cubed)
{1: 1, -1/2 - sqrt(3)*I/2: 1, -1/2 + sqrt(3)*I/2: 1}
>>> # real_roots() returns only real roots
>>> real_roots(cubed)
[1]
real_roots() 调用 RootOf()，因此对于所有根为实数的方程，通过迭代方程的根数，可以得到相同的结果：
>>> [RootOf(expression, n) for n in range(3)]
[-2, -2, 3]
近似数值解与根重数
nroots
nroots() 给出多项式根的近似数值解。此示例表明它可能包含数值噪声，例如本应是实根的部分（可忽略的）虚部分：
>>> nroots(expression)
[3.0, -2.0 - 4.18482169793536e-14*I, -2.0 + 4.55872552179222e-14*I]
如果你想要实根的数值近似，但又想确切知道哪些根是实数，那么最好的方法是使用 real_roots() 结合 evalf()：
>>> [r.n(2) for r in real_roots(expression)]
[-2.0, -2.0, 3.0]
>>> [r.is_real for r in real_roots(expression)]
[True, True, True]
nroots() 类似于 NumPy 的 roots() 函数。通常这两者的区别在于 nroots() 更精确但速度较慢。
nroots() 的一个主要优势是它可以计算任何多项式的数值近似根，只要其系数可以通过 evalf() 进行数值评估（即它们没有自由符号）。相反，符号解可能对于高阶（五阶或更高阶）多项式不可能，正如 阿贝尔-鲁菲尼定理 所解释的那样。即使存在闭合形式解，它们可能有太多项以至于在实际中不太有用。因此，即使存在闭合形式的符号解，你可能还是会选择使用 nroots() 来找到近似的数值解。例如，四阶（四次）多项式的闭合形式根可能会相当复杂：
>>> rq0, rq1, rq2, rq3 = roots(x**4 + 3*x**2 + 2*x + 1)
>>> rq0
sqrt(-4 - 2*(-1/8 + sqrt(237)*I/36)**(1/3) + 4/sqrt(-2 + 7/(6*(-1/8 + sqrt(237)*I/36)**(1/3)) + 2*(-1/8 + sqrt(237)*I/36)**(1/3)) - 7/(6*(-1/8 + sqrt(237)*I/36)**(1/3)))/2 - sqrt(-2 + 7/(6*(-1/8 + sqrt(237)*I/36)**(1/3)) + 2*(-1/8 + sqrt(237)*I/36)**(1/3))/2
因此，你可能更喜欢一个近似数值解：
>>> rq0.n()
-0.349745826211722 - 0.438990337475312*I
nroots() 有时对于数值条件不佳的多项式可能会失败，例如 威尔金森多项式。使用 RootOf() 和 evalf()，如 数值评估 CRootOf 的根 中描述的，可以避免由于使用更精确但速度较慢的数值算法（基于孤立区间）而导致的病态和返回虚假复数部分。
复数根
对于复数根，可以使用类似的函数，例如 solve()：
>>> from sympy import solve, roots, nroots, real_roots, expand, RootOf, CRootOf, Symbol
>>> from sympy import Poly
>>> from sympy.abc import x
>>> expression_complex = (x**2+4)**2 * (x-3)
>>> solve(expression_complex, x, dict=True)
[{x: 3}, {x: -2*I}, {x: 2*I}]
如果常数是符号性的，你可能需要指定它们的域以便于 SymPy 认识到解不是实数。例如，指定 (a) 为正会导致虚数根：
>>> a = Symbol("a", positive=True)
>>> symbolic_complex = (x**2+a)**2 * (x-3)
>>> solve(symbolic_complex, x, dict=True)
[{x: 3}, {x: -I*sqrt(a)}, {x: I*sqrt(a)}]
roots() 也可以找到虚根或复根：
>>> roots(expression_complex, x)
{3: 1, -2*I: 2, 2*I: 2}
RootOf() 也会返回复杂根：
>>> [RootOf(expression_complex, n) for n in range(0,3)]
[3, -2*I, -2*I]
real_roots() 将仅返回实根。
>>> real_roots(expression_complex)
[3]
real_roots() 的优点在于，它可能比生成所有根更有效：RootOf() 对于复杂根可能会比较慢。
如果您将表达式转换为多项式类 Poly，则可以使用其 all_roots() 方法查找根：
>>> expression_complex_poly = Poly(expression_complex)
>>> expression_complex_poly.all_roots()
[3, -2*I, -2*I, 2*I, 2*I]
使用解决方案结果
从结果中提取解的方式取决于结果的形式。
列表（all_roots, real_roots, nroots）
您可以使用标准的 Python 列表遍历技术进行遍历。在这里，我们将每个根代入表达式中以验证结果为 (0)：
>>> expression = (x+2)**2 * (x-3)
>>> my_real_roots = real_roots(expression)
>>> my_real_roots
[-2, -2, 3]
>>> for root in my_real_roots:
... print(f"expression({root}) = {expression.subs(x, root)}")
expression(-2) = 0
expression(-2) = 0
expression(3) = 0
字典列表（solve）
请参考 使用解决方案结果。
字典（roots）
您可以使用标准的 Python 列表遍历技术，比如遍历字典中的键和值。这里我们打印每个根的值和重复次数：
>>> my_roots = roots(expression)
>>> my_roots
{-2: 2, 3: 1}
>>> for root, multiplicity in my_roots.items():
... print(f"Root {root} has multiplicity of {multiplicity}")
Root 3 has multiplicity of 1
Root -2 has multiplicity of 2
表达式（factor）
您可以使用各种 SymPy 技术来操作代数表达式，例如用符号或数值替换 (x)：
>>> from sympy.abc import y
>>> factored = factor(expression_expanded)
>>> factored
(x - 3)*(x + 2)**2
>>> factored.subs(x, 2*y)
(2*y - 3)*(2*y + 2)**2
>>> factored.subs(x, 7)
324
折衷方案
数学精确性、根列表的完整性和速度
考虑高阶多项式 (x⁵ - x + 1 = 0)。nroots() 返回所有五个根的数值近似：
>>> from sympy import roots, solve, real_roots, nroots
>>> from sympy.abc import x
>>> fifth_order = x**5 - x + 1
>>> nroots(fifth_order)
[-1.16730397826142,
-0.181232444469875 - 1.08395410131771*I,
-0.181232444469875 + 1.08395410131771*I,
0.764884433600585 - 0.352471546031726*I,
0.764884433600585 + 0.352471546031726*I]
roots() 有时可能只返回部分根，或者如果无法用根式表示任何根，则返回空集合：
>>> roots(fifth_order, x)
{}
但如果您设置标志 strict=True，roots() 将告知您无法返回所有根：
>>> roots(x**5 - x + 1, x, strict=True)
Traceback (most recent call last):
...
sympy.polys.polyerrors.UnsolvableFactorError: Strict mode: some factors cannot be solved in radicals, so a complete
list of solutions cannot be returned. Call roots with strict=False to
get solutions expressible in radicals (if there are any).
获取所有根，也许是隐含的
solve() 将作为 CRootOf（ComplexRootOf()）类成员返回所有五个根
>>> fifth_order_solved = solve(fifth_order, x, dict=True)
>>> fifth_order_solved
[{x: CRootOf(x**5 - x + 1, 0)},
{x: CRootOf(x**5 - x + 1, 1)},
{x: CRootOf(x**5 - x + 1, 2)},
{x: CRootOf(x**5 - x + 1, 3)},
{x: CRootOf(x**5 - x + 1, 4)}]
每个 CRootOf 中的第二个参数是根的索引。
数值评估 CRootOf 根
您可以使用 evalf() 对那些 CRootOf 根进行数值评估：
>>> for root in fifth_order_solved:
... print(root[x].n(10))
-1.167303978
-0.1812324445 - 1.083954101*I
-0.1812324445 + 1.083954101*I
0.7648844336 - 0.352471546*I
0.7648844336 + 0.352471546*I
如果您只对唯一的实根感兴趣，使用 real_roots() 更快，因为它不会尝试找到复数根：
>>> real_root = real_roots(fifth_order, x)
>>> real_root
[CRootOf(x**5 - x + 1, 0)]
>>> real_root[0].n(10)
-1.167303978
表示根
RootOf()、real_roots() 和 all_roots() 可以精确地找到多项式的所有根，尽管存在 Abel-Ruffini 定理。这些函数允许精确地分类和符号操作根。
>>> from sympy import init_printing
>>> init_printing()
>>> real_roots(fifth_order)
/ 5 \
[CRootOf\x - x + 1, 0/]
>>> r = r0, r1, r2, r3, r4 = Poly(fifth_order, x).all_roots(); r
/ 5 \ / 5 \ / 5 \ / 5 \ / 5 \
[CRootOf\x - x + 1, 0/, CRootOf\x - x + 1, 1/, CRootOf\x - x + 1, 2/, CRootOf\x - x + 1, 3/, CRootOf\x - x + 1, 4/]
>>> r0
/ 5 \
CRootOf\x - x + 1, 0/
现在已经精确找到了根，可以确定它们的属性而不受数值噪声的影响。例如，我们可以判断根是实数还是虚数。例如，如果我们请求根 r1 的 conjugate()（实部相同，虚部相反），并且这恰好等于另一个根 r2，则根 r2 将被返回：
>>> r0.n()
-1.16730397826142
>>> r0.is_real
True
>>> r1.n()
-0.181232444469875 - 1.08395410131771*I
>>> r2.n()
-0.181232444469875 + 1.08395410131771*I
>>> r1
/ 5 \
CRootOf\x - x + 1, 1/
>>> r1.conjugate()
/ 5 \
CRootOf\x - x + 1, 2/
>>> r1.is_real
False
solve() 在可能的情况下也会给出复数根，但比直接使用 all_roots() 效率低。
RootOf() 以可以符号操作和任意精度计算的方式精确表示根。RootOf() 的表示使得能够精确地：
-
计算具有精确有理系数的多项式的所有根。
-
确定每个根的多重性。
-
精确确定根是否为实数。
-
精确地排序实根和复根。
-
知道哪些根是复共轭对。
-
确切确定哪些根是有理数，哪些是无理数。
-
精确地表示每个可能的代数数。
其他数值方法如 NumPy 的roots()，nroots()和nsolve()在所有情况下都无法稳健地执行这些操作。同样地，当使用evalf()进行数值评估时，由solve()或roots()返回的根式表达式也无法稳健地执行这些操作。
并非所有方程都能求解
没有闭合形式解的方程
如上所述，高阶多项式（五次或更高次）不太可能有闭合形式的解，因此你可能需要用例如RootOf 如上所述，或使用数值方法如nroots 如上所述来表示它们。
报告错误
如果您在使用这些命令时遇到错误，请在SymPy 邮件列表上发布问题。在问题得到解决之前，您可以使用另一个寻找多项式根的函数或尝试考虑的替代方案之一。
代数求解一个矩阵方程
原文：
docs.sympy.org/latest/guides/solving/solve-matrix-equation.html
使用 SymPy 解矩阵（线性）方程。例如，解 ( \left[\begin{array}{cc} c & d\1 & -e\end{array}\right] \left[\begin{array}{cc} x\y\end{array}\right] = \left[\begin{array}{cc} 2\0\end{array}\right] ) 得到 ( \left[\begin{array}{cc} x\y\end{array}\right] = \left[\begin{array}{cc} \frac{2e}{ce+d}\\frac{2}{ce+d}\end{array}\right])。
可供考虑的替代方法
-
如果你的矩阵和常数向量只包含数字，而不是符号，例如 (\left[\begin{array}{cc} 1 & 2\3 & 4\end{array}\right] \left[\begin{array}{cc} x\y\end{array}\right] = \left[\begin{array}{cc} 2\0\end{array}\right])，你可以使用 SymPy 的其他免费开源软件包之一，而不是 SymPy：
-
NumPy 的
numpy.linalg.solve() -
SciPy 的
scipy.linalg.solve() -
mpmath 的 lu_solve()
-
-
解矩阵方程等同于解线性方程组，所以如果你愿意，你可以代数求解一个线性方程组
-
如果你已经将问题表述为一个线性方程组，并希望将其转换为矩阵形式，可以使用
linear_eq_to_matrix()函数，然后按照本指南的步骤进行操作。
求解矩阵方程
这里是使用 SymPy 的 sympy.matrices.matrixbase.MatrixBase.solve() 求解矩阵方程的示例。我们使用标准的矩阵方程形式 (Ax=b)，其中
-
(A) 是表示线性方程中系数的矩阵
-
(x) 是要求解的未知数的列向量
-
(b) 是常数的列向量，其中每行是一个方程的值
>>> from sympy import init_printing
>>> init_printing(use_unicode=True)
>>> from sympy import symbols
>>> from sympy.matrices import Matrix
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c,d], [1, -e]])
>>> A
⎡c d ⎤
⎢ ⎥
⎣1 -e⎦
>>> b = Matrix([2, 0])
>>> b
⎡2⎤
⎢ ⎥
⎣0⎦
>>> A.solve(b)
⎡ 2⋅e ⎤
⎢───────⎥
⎢c⋅e + d⎥
⎢ ⎥
⎢ 2 ⎥
⎢───────⎥
⎣c⋅e + d⎦
指南
矩阵通常必须是方阵
矩阵 (A) 通常必须是方阵，以表示具有与方程数量相同的未知数的线性方程组。如果不是，SymPy 将会报错ShapeError: `self` and `rhs` must have the same number of rows.
例外的是，SymPy 使用Moore-Penrose 伪逆 的要求，这不需要矩阵是方阵。
解矩阵方程的方法
SymPy 的矩阵求解方法，sympy.matrices.matrixbase.MatrixBase.solve()，可以使用几种不同的方法，这些方法在 API 参考链接中列出。根据矩阵的性质，某种方法可能更有效。默认情况下，将使用高斯-约当消元法。
在 solve 中指定一个方法相当于使用专门的求解函数。例如，使用solve和method='LU'调用LUsolve()。
解决多个相同矩阵方程
如果需要重复解决具有相同矩阵(A)但不同常向量(b)的矩阵方程，则更有效的方法是使用以下方法之一。
您可以通过LUsolve()使用LU 分解：
>>> from sympy import symbols, Matrix, eye, simplify
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c,d], [1, -e]])
>>> A
⎡c d ⎤
⎢ ⎥
⎣1 -e⎦
>>> b = Matrix([2, 0])
>>> b
⎡2⎤
⎢ ⎥
⎣0⎦
>>> solution = A.LUsolve(b)
>>> solution
⎡ 2⋅e ⎤
⎢───────⎥
⎢c⋅e + d⎥
⎢ ⎥
⎢ 2 ⎥
⎢───────⎥
⎣c⋅e + d⎦
>>> # Demonstrate that solution is correct
>>> simplify(A * solution)
⎡2⎤
⎢ ⎥
⎣0⎦
>>> b2 = Matrix([4, 0])
>>> b2
⎡4⎤
⎢ ⎥
⎣0⎦
>>> solution2 = A.LUsolve(b2)
>>> solution2
⎡ 4⋅e ⎤
⎢───────⎥
⎢c⋅e + d⎥
⎢ ⎥
⎢ 4 ⎥
⎢───────⎥
⎣c⋅e + d⎦
>>> # Demonstrate that solution2 is correct
>>> simplify(A * solution2)
⎡4⎤
⎢ ⎥
⎣0⎦
另一种方法是计算逆矩阵，但这几乎总是比较慢，对于更大的矩阵来说，速度慢得多。如果高效计算不是优先考虑的话，可以使用inv():
>>> from sympy import symbols, Matrix, simplify
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c,d], [1, -e]])
>>> b = Matrix([2, 0])
>>> b
⎡2⎤
⎢ ⎥
⎣0⎦
>>> b2 = Matrix([4, 0])
>>> b2
⎡4⎤
⎢ ⎥
⎣0⎦
>>> inv = A.inv()
>>> inv
⎡ e d ⎤
⎢─────── ───────⎥
⎢c⋅e + d c⋅e + d⎥
⎢ ⎥
⎢ 1 -c ⎥
⎢─────── ───────⎥
⎣c⋅e + d c⋅e + d⎦
>>> # Solves Ax = b for x
>>> solution = inv * b
>>> solution
⎡ 2⋅e ⎤
⎢───────⎥
⎢c⋅e + d⎥
⎢ ⎥
⎢ 2 ⎥
⎢───────⎥
⎣c⋅e + d⎦
>>> # Demonstrate that solution is correct
>>> simplify(A * solution)
⎡2⎤
⎢ ⎥
⎣0⎦
>>> # Solves Ax = b2 for x
>>> solution2 = inv * b2
>>> solution2
⎡ 4⋅e ⎤
⎢───────⎥
⎢c⋅e + d⎥
⎢ ⎥
⎢ 4 ⎥
⎢───────⎥
⎣c⋅e + d⎦
>>> # Demonstrate that solution2 is correct
>>> simplify(A * solution2)
⎡4⎤
⎢ ⎥
⎣0⎦
确定大型符号矩阵的逆可能不可计算。
使用符号矩阵
操作符号矩阵的计算复杂性随着矩阵大小的增加而迅速增加。例如，符号矩阵行列式中的项数随矩阵维数的阶乘增加。因此，可以解决的矩阵的最大维度比数值矩阵更有限。例如，这个 4x4 符号矩阵的行列式有 24 个项，每个项有四个元素：
>>> from sympy import MatrixSymbol
>>> A = MatrixSymbol('A', 4, 4).as_explicit()
>>> A
⎡A₀₀ A₀₁ A₀₂ A₀₃⎤
⎢ ⎥
⎢A₁₀ A₁₁ A₁₂ A₁₃⎥
⎢ ⎥
⎢A₂₀ A₂₁ A₂₂ A₂₃⎥
⎢ ⎥
⎣A₃₀ A₃₁ A₃₂ A₃₃⎦
>>> A.det()
A₀₀⋅A₁₁⋅A₂₂⋅A₃₃ - A₀₀⋅A₁₁⋅A₂₃⋅A₃₂ - A₀₀⋅A₁₂⋅A₂₁⋅A₃₃ + A₀₀⋅A₁₂⋅A₂₃⋅A₃₁ +
A₀₀⋅A₁₃⋅A₂₁⋅A₃₂ - A₀₀⋅A₁₃⋅A₂₂⋅A₃₁ - A₀₁⋅A₁₀⋅A₂₂⋅A₃₃ + A₀₁⋅A₁₀⋅A₂₃⋅A₃₂ +
A₀₁⋅A₁₂⋅A₂₀⋅A₃₃ - A₀₁⋅A₁₂⋅A₂₃⋅A₃₀ - A₀₁⋅A₁₃⋅A₂₀⋅A₃₂ + A₀₁⋅A₁₃⋅A₂₂⋅A₃₀ +
A₀₂⋅A₁₀⋅A₂₁⋅A₃₃ - A₀₂⋅A₁₀⋅A₂₃⋅A₃₁ - A₀₂⋅A₁₁⋅A₂₀⋅A₃₃ + A₀₂⋅A₁₁⋅A₂₃⋅A₃₀ +
A₀₂⋅A₁₃⋅A₂₀⋅A₃₁ - A₀₂⋅A₁₃⋅A₂₁⋅A₃₀ - A₀₃⋅A₁₀⋅A₂₁⋅A₃₂ + A₀₃⋅A₁₀⋅A₂₂⋅A₃₁ +
A₀₃⋅A₁₁⋅A₂₀⋅A₃₂ - A₀₃⋅A₁₁⋅A₂₂⋅A₃₀ - A₀₃⋅A₁₂⋅A₂₀⋅A₃₁ + A₀₃⋅A₁₂⋅A₂₁⋅A₃₀
并解决它大约需要一分钟，而类似的 3x3 矩阵则少于一秒。矩阵中不相关的符号条目越多，操作速度就越慢。这个例子是在所有元素都是独立符号的矩阵中找到一个通用解的极端情况，因此对于其大小而言是最慢的。
加速解决矩阵方程
这里有一些建议：
-
如果矩阵元素为零，请确保它们被识别为零。您可以通过将它们设为零或应用假设来实现这一点。
-
选择适合于矩阵性质的求解方法，例如埃尔米特、对称或三角形式。参见解矩阵方程的方法。
-
使用
DomainMatrix类，它可能更快，因为它限制了矩阵元素的定义域。
使用解的结果
将解作为向量使用
您可以将解结果用作向量。例如，为了证明解 (x) 是正确的，您可以将其与矩阵 (A) 相乘，并验证其是否生成常数向量 (b)：
>>> from sympy import symbols, simplify
>>> from sympy.matrices import Matrix
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c,d], [1, -e]])
>>> b = Matrix([2, 0])
>>> solution = A.solve(b)
>>> solution
⎡ 2⋅e ⎤
⎢───────⎥
⎢c⋅e + d⎥
⎢ ⎥
⎢ 2 ⎥
⎢───────⎥
⎣c⋅e + d⎦
>>> # Not immediately obvious whether this result is a zeroes vector
>>> (A * solution) - b
⎡ 2⋅c⋅e 2⋅d ⎤
⎢─────── + ─────── - 2⎥
⎢c⋅e + d c⋅e + d ⎥
⎢ ⎥
⎣ 0 ⎦
>>> # simplify reveals that this result is a zeroes vector
>>> simplify((A * solution) - b)
⎡0⎤
⎢ ⎥
⎣0⎦
请注意，我们不得不使用 simplify() 来使 SymPy 简化矩阵元素中的表达式，以便立即明确解是正确的。
从解中提取元素
因为您可以通过遍历列向量中的元素，可以使用标准的 Python 技术提取其元素。例如，您可以使用列表推导创建元素列表
>>> [element for element in solution]
⎡ 2⋅e 2 ⎤
⎢───────, ───────⎥
⎣c⋅e + d c⋅e + d⎦
或者可以通过下标提取单个元素
>>> solution[0]
2⋅e
───────
c⋅e + d
方程无解
如果矩阵的行列式为零，则与之相关的矩阵方程无解：
>>> from sympy import symbols
>>> from sympy.matrices import Matrix
>>> c, d, e = symbols("c, d, e")
>>> A = Matrix([[c*e**2, d*e], [c*e, d]])
>>> A
⎡ 2 ⎤
⎢c⋅e d⋅e⎥
⎢ ⎥
⎣c⋅e d ⎦
>>> b = Matrix([2, 0])
>>> A.LUsolve(b)
Traceback (most recent call last):
...
NonInvertibleMatrixError: Matrix det == 0; not invertible.
报告错误
如果您在矩阵求解函数中发现错误，请在SymPy 邮件列表上发布问题。在问题解决之前，您可以考虑使用考虑的替代方法中列出的其他方法。
浙公网安备 33010602011771号