SymPy 1.13 中文文档(五)
代数方式解方程
原文:
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 邮件列表上发布问题。在问题解决之前,您可以考虑使用考虑的替代方法中列出的其他方法。
单变量代数系统的一个或一组不等式简化
原文:
docs.sympy.org/latest/guides/solving/reduce-inequalities-algebraically.html
使用 SymPy 在单变量代数中简化一个或一组不等式。例如,简化 (x² < \pi),(x > 0) 将得到 (0 < x < \sqrt{\pi})。
注意
SymPy 目前仅能简化不等式中的一个符号(变量)。
SymPy 可以简化包含多个符号的系统,如果每个不等式只有一个符号。
考虑的替代方案
-
若要简化不等式中的多个符号,请尝试使用 SciPy 的
linprog()
-
要简化布尔表达式,请使用
as_set
例子
简化单变量代数不等式系统
reduce_inequalities()
接受要作为系统简化的不等式列表或元组:
>>> from sympy import symbols, reduce_inequalities, pi
>>> x = symbols('x')
>>> reduce_inequalities([x >= 0, x**2 <= pi], x)
(0 <= x) & (x <= sqrt(pi))
注意
虽然solve()
目前可以通过在内部调用reduce_inequalities()
来完成相同的功能,但该功能可能会在solve()
中被弃用或删除。因此,我们建议使用reduce_inequalities()
。
reduce_inequalities()
是顶层不等式简化函数,将在需要时内部调用任何其他低级不等式简化函数。
简化单变量代数不等式系统
如果只有一个不等式,可以选择排除列表结构,并将reduce_inequalities()
作为表达式传递给它:
>>> from sympy import symbols, reduce_inequalities, pi
>>> x = symbols('x')
>>> reduce_inequalities(x**2 <= pi, x)
(x <= sqrt(pi)) & (-sqrt(pi) <= x)
指南
在函数调用中包含要简化的变量
我们建议您将要简化的变量作为reduce_inequalities()
的第二个参数,以确保它对所需变量进行简化。
代数地减少一组不等式
您可以创建您的不等式,然后将系统简化为列表:
>>> from sympy import symbols, reduce_inequalities, pi
>>> x = symbols('x')
>>> reduce_inequalities([3*x >= 1, x**2 <= pi], x)
(1/3 <= x) & (x <= sqrt(pi))
使用结果
使用结果的常见方式是提取符号(变量)的边界。例如,对于 (0 < x < \sqrt{\pi}) 的解,您可能希望提取 (0) 和 (\sqrt{\pi})。
提取分解关系列表
您可以将通过 ^
(Or
) 或 &
(And
) 连接的一组关系分解为单个关系使用关系原子。使用canonical
将为每个关系放置顺序,使符号在左侧,因此您可以获取右侧rhs
以提取常数:
>>> from sympy import symbols, reduce_inequalities, pi
>>> from sympy.core.relational import Relational
>>> x = symbols('x')
>>> eq = reduce_inequalities([3*x >= 1, x**2 <= pi], x); eq
(1/3 <= x) & (x <= sqrt(pi))
>>> relations = [(i.lhs, i.rel_op, i.rhs) for i in [i.canonical for i in eq.atoms(Relational)]]
>>> relations_sorted = sorted(relations, key=lambda x: float(x[2])) # Sorting relations just to ensure consistent list order for docstring testing
>>> relations_sorted
[(x, '>=', 1/3), (x, '<=', sqrt(pi))]
提取关系元组
简化关系的args
(参数)是单独的关系,因此您可以从左侧或右侧的args
中提取常数:
>>> from sympy import symbols, reduce_inequalities, pi
>>> x = symbols('x')
>>> eq = reduce_inequalities([3*x >= 1, x**2 <= pi], x); eq
(1/3 <= x) & (x <= sqrt(pi))
>>> eq.args
(1/3 <= x, x <= sqrt(pi))
>>> constants = []
>>> for arg in eq.args:
... if arg.lhs == x:
... constants.append(arg.rhs)
... else:
... constants.append(arg.lhs)
>>> constants
[1/3, sqrt(pi)]
使用 SymPy 减少不等式的限制
SymPy 只能对感兴趣的每个不等式中的一个符号进行简化。
SymPy 目前只能针对给定不等式中感兴趣的一个符号(变量)进行简化。
>>> from sympy import reduce_inequalities, symbols
>>> x, y = symbols("x y")
>>> reduce_inequalities([x + y > 1, y > 0], [x, y])
Traceback (most recent call last):
...
NotImplementedError: inequality has more than one symbol of interest.
使用 SciPy 的linprog()
可以减少这个不等式系统。
SymPy 可以在系统中对超过一个符号进行简化,如果每个不等式只有一个感兴趣的符号。例如,以下不等式系统包含两个变量,(x) 和 (y)。SymPy 可以对 (x) 进行简化,并给出 (y) 的约束条件。
>>> from sympy import reduce_inequalities, symbols
>>> x, y = symbols("x y")
>>> reduce_inequalities([x + y > 1, y > 0], x)
(0 < y) & (y < oo) & (x > 1 - y)
(oo
是Infinity
.)
如果每个不等式仅包含一个要简化的符号,SymPy 可以为多个符号减少不等式集合:
>>> from sympy import reduce_inequalities, symbols
>>> x, y = symbols("x y")
>>> x_y_reduced = reduce_inequalities([x > 1, y > 0], [x, y]); x_y_reduced
(0 < y) & (1 < x) & (x < oo) & (y < oo)
请注意,这提供的数学洞察力仅限于分别减少不等式:
>>> from sympy import And
>>> x_reduced = reduce_inequalities(x > 1, x); x_reduced
(1 < x) & (x < oo)
>>> y_reduced = reduce_inequalities(y > 0, y); y_reduced
(0 < y) & (y < oo)
>>> And(x_reduced, y_reduced) == x_y_reduced
True
因此,解决此类不等式作为集合的好处可能只是方便性。
SymPy 能够解决的不等式类型限制
reduce_inequalities()
可以解决涉及要简化符号的幂或涉及另一个符号的不等式系统:
>>> from sympy import reduce_inequalities
>>> from sympy.abc import x, y
>>> reduce_inequalities([x ** 2 < 4, x > 0], x)
(0 < x) & (x < 2)
>>> reduce_inequalities([x < y, x > 0], x)
(0 < x) & (x < oo) & (x < y)
>>> reduce_inequalities([x ** 2 - y < 4, x > 0], x)
Traceback (most recent call last):
...
NotImplementedError: The inequality, -_y + x**2 - 4 < 0, cannot be solved using
solve_univariate_inequality.
并非所有周期函数的结果都会被返回
对于三角不等式返回的结果受其周期间隔的限制。reduce_inequalities()
试图返回足够的解,以便所有(无限多个)解都可以通过返回的解加上方程的整数倍的 periodicity()
(这里是 (2\pi))生成。
>>> from sympy import reduce_inequalities, cos
>>> from sympy.abc import x, y
>>> from sympy.calculus.util import periodicity
>>> reduce_inequalities([2*cos(x) < 1, x > 0], x)
(0 < x) & (x < oo) & (pi/3 < x) & (x < 5*pi/3)
>>> periodicity(2*cos(x), x)
2*pi
并非所有不等式系统都可以简化
无法满足的不等式系统
如果不等式系统具有不兼容的条件,例如 (x < 0) 和 (x > \pi),SymPy 将返回 False
:
>>> from sympy import symbols, reduce_inequalities, pi
>>> x = symbols('x')
>>> reduce_inequalities([x < 0, x > pi], x)
False
无法在解析上简化的不等式系统
SymPy 可能会反映您的不等式系统在代数(符号)上无法表达的解不存在,如返回诸如 NotImplementedError
的错误:
>>> from sympy import symbols, reduce_inequalities, cos
>>> x = symbols('x')
>>> reduce_inequalities([cos(x) - x > 0, x > 0], x)
Traceback (most recent call last):
...
NotImplementedError: The inequality, -x + cos(x) > 0, cannot be solved using solve_univariate_inequality.
因此,您可能需要使用 SciPy 的 linprog()
在数值上简化您的不等式。
可以在解析上简化的不等式,但 SymPy 无法简化的系统
请参阅上文的 使用 SymPy 进行不等式简化的限制。
报告 Bug
如果您在 diophantine()
中发现 Bug,请在 SymPy 邮件列表 上发布问题。在问题解决之前,您可以使用列在 考虑的替代方法 中的其他方法。