Sympy "wrongly" 求解 ODE 的系统

Question

我在 sympy 中定义了两个 ODE，如下所示：

import sympy as sp

a = sp.Symbol("a", positive=True)
b = sp.Symbol("b", positive=True)
t = sp.Symbol("t")
x = sp.Function("x")
y = sp.Function("y")

eq = (
    sp.Eq(sp.Derivative(x(t), t), -a * x(t)),
    sp.Eq(sp.Derivative(y(t), t), a * x(t) - b * y(t)),
)

当我尝试使用 dsolve 解决它们时，我得到以下结果：

>>> sp.dsolve(eq)

[Eq(x(t), -C1*(a - b)*exp(-a*t)/a),
 Eq(y(t), C1*exp(-a*t) + C2*exp(-b*t))]

单独求解第一个方程时，我得到以下结果：

>>> sp.dsolve(eq[0])

Eq(x(t), C1*exp(-a*t))

y(t) 中不依赖于第二个 ODE 的第一个 ODE 怎么可能突然包含 b？用 mathematica 解决它似乎给出了正确的结果：

x(t) = C1*exp(-at)
y(t) = C2*exp(-bt) - (a*exp(-at-bt)*(-exp(at)+exp(bt))*C1)/(a-b)

任何帮助 python 解决此问题的人都将不胜感激。

Answer 1

SymPy 的解决方案是正确的。它也等同于 Mathematica 中的一个，只是任意常数的定义不同。您可以替换 C1 -> -C1*(a/(a-b)) 之类的内容以从一种形式转换为另一种形式：

In [49]: sp.dsolve(eq)
Out[49]: 
⎡                    -a⋅t                             ⎤
⎢       -C₁⋅(a - b)⋅ℯ                  -a⋅t       -b⋅t⎥
⎢x(t) = ──────────────────, y(t) = C₁⋅ℯ     + C₂⋅ℯ    ⎥
⎣               a                                     ⎦

如果您改变要求解的函数的顺序，您还可以获得稍微不同的形式：

In [53]: sp.dsolve(eq, [y(t), x(t)])
Out[53]: 
⎡               -a⋅t                            ⎤
⎢         C₁⋅a⋅ℯ           -b⋅t             -a⋅t⎥
⎢y(t) = - ────────── + C₂⋅ℯ    , x(t) = C₁⋅ℯ    ⎥
⎣           a - b                               ⎦

在每个方程中单独exp(-a*t)的系数可以是任意常数。通常其他参数如 a 和 b 可以吸收到任意常数中以简化结果。但是，x 和 y 的解中的系数之间存在约束，因此一个方程的系数必须与另一个方程中的系数相差 -a/(a-b)。这里 SymPy 选择吸收最后一个方程中的参数，但这意味着它们必须重新出现在第一个方程中。

您可以在 Mathematica 显示的输出中看到相同的因子：

⎡                               ⎛   a⋅t    b⋅t⎞  -a⋅t - b⋅t           ⎤
⎢           -a⋅t           C₁⋅a⋅⎝- ℯ    + ℯ   ⎠⋅ℯ                 -b⋅t⎥
⎢x(t) = C₁⋅ℯ    , y(t) = - ──────────────────────────────── + C₂⋅ℯ    ⎥
⎣                                       a - b                         ⎦

所有这些都是原始 ODE 的有效通用解决方案，但在任意常数的定义方式上有所不同，因此这只是您喜欢的口味问题。

我个人最不喜欢 Mathematica 输出，因为它混合了两个指数项的系数。我可以理解它的来源，因为它基本上是矩阵指数的原始输出：

In [60]: M = sp.Matrix([[-a, 0], [a, -b]])

In [60]: (M*t).exp().applyfunc(sp.factor) * sp.Matrix([C1, C2])
Out[60]: 
⎡                    -a⋅t                 ⎤
⎢                C₁⋅ℯ                     ⎥
⎢                                         ⎥
⎢     ⎛ a⋅t    b⋅t⎞  -a⋅t  -b⋅t           ⎥
⎢C₁⋅a⋅⎝ℯ    - ℯ   ⎠⋅ℯ    ⋅ℯ           -b⋅t⎥
⎢────────────────────────────── + C₂⋅ℯ    ⎥
⎣            a - b                        ⎦

SymPy 公式略有不同，使用 Jordan 形式的矩阵指数隐式吸收特征向量矩阵 inv(P) 为常数：

In [63]: P, J = M.jordan_form()

In [64]: P*(J*t).exp() * sp.Matrix([C1, C2])
Out[64]: 
⎡    ⎛     b⎞  -a⋅t ⎤
⎢ C₁⋅⎜-1 + ─⎟⋅ℯ     ⎥
⎢    ⎝     a⎠       ⎥
⎢                   ⎥
⎢    -a⋅t       -b⋅t⎥
⎣C₁⋅ℯ     + C₂⋅ℯ    ⎦

Sympy "wrongly" 求解 ODE 的系统

Sympy "wrongly" solving system of ODEs

python

sympy