代数求解常微分方程 (ODE)¶
使用SymPy代数求解常微分方程（ODE）。例如，求解 \(y''(x) + 9y(x)=0 \) 得到 \( y(x)=C_{1} \sin(3x)+ C_{2} \cos(3x)\)。
考虑的替代方案¶
要数值求解一个常微分方程组，可以使用 SciPy ODE 求解器 如
solve_ivp。你也可以使用 SymPy 创建一个常微分方程，然后lambdify()它，以便使用 SciPy 的solve_ivp进行数值求解，如下文 Numerically Solve an ODE in SciPy 所述。
求解常微分方程 (ODE)¶
以下是使用 dsolve() 代数求解上述常微分方程的示例。然后，您可以使用 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)
toctree 是一个 reStructuredText :dfn:指令 ，这是一个非常多功能的标记。指令可以有参数、选项和内容。
指导¶
定义导数¶
表达函数导数的方法有很多。对于一个未定义的函数，Derivative 和 diff() 都表示未定义的导数。因此，以下所有的 ypp（“y 双撇”）都表示 \(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} 等形式给出，例如，\(f(x)\) 在 \(x = x_{0}\) 处的值为 \(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\) 已经将其独立变量作为参数 \(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)
使用任意常量¶
你可以操作诸如 C1、C2 和 C3 这样的任意常量，这些常量是由 dsolve() 自动生成的，方法是将它们创建为符号。例如，如果你想为任意常量赋值，你可以将它们创建为符号，然后使用 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 中数值求解常微分方程¶
利用 SciPy 快速数值求解常微分方程的常见工作流程是
在 SymPy 中设置一个 ODE
使用
lambdify()将其转换为数值函数通过 使用 SciPy 的
solve_ivp对 ODE 进行数值积分 来解决初值问题。
以下是一个来自化学动力学领域的例子，其中非线性常微分方程采用这种形式：
和
>>> 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()
SciPy 的 solve_ivp 返回一个结果，其中包含 y（数值函数结果，这里为浓度）值，对应于每个化学物种，以及对应的时间点 t_eval。
常微分方程求解提示¶
返回未计算的积分¶
默认情况下，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))
无闭式解的方程¶
如上所述， Some ODEs Cannot Be Solved Explicitly, Only Implicitly。
此外，一些微分方程系统没有封闭形式的解，因为它们是混沌的，例如 Lorenz 系统 或由这两个微分方程描述的双摆（简化为 ScienceWorld）：
>>> 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
对于这种情况，你可以按照 Alternatives to Consider 中提到的方法，通过数值方法求解方程。
报告一个错误¶
如果你在使用 dsolve() 时发现了一个错误，请在 SymPy 邮件列表 上发布这个问题。在问题解决之前，你可以使用 Alternatives to Consider 中列出的其他方法。