SymPy 1.13 中文文档(四)
编写自定义函数
本指南将描述如何在 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^{-2t} & \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)作为装置(参考下面的框图)。
解决方案
>>> # 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 邮件列表上发布问题。