Sympy使用dsolve求解常微分方程组结果异常问题咨询

Sympy求解ODE方程组时使用全局统一的任意常数，和单方程求解时的局部常数定义不同，你看到的x(t)表达式含b并非求解错误，只是常数未做归一化导致的形式差异：

• 联立求解返回的x(t) = -C1*(a - b)*exp(-a*t)/a中，-C1*(a-b)/a整体就是一个不随t变化的任意常数，完全可以等价替换为单方程求解时的常数C1，b只是常数系数的一部分，不会改变x(t)按exp(-a*t)衰减的演化规律。
• 对比Mathematica的输出，两者的y(t)本质也是等价的：把Mathematica给出的y(t)展开化简后，和Sympy返回的y(t) = C1*exp(-a*t) + C2*exp(-b*t)仅差常数的线性组合，数学上完全等价。

可以直接用sp.checkodesol(eq, sol)验证Sympy原始返回的解，返回结果为True就说明解满足原方程，不存在计算错误。

可以用两种方式得到和单方程求解、Mathematica输出形式一致的结果：

• 顺序求解法：先解独立的x(t)方程，再代入y(t)的方程求解，从根源避免全局常数带来的系数冗余

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)),
)
# 先求解独立的x方程
x_sol = sp.dsolve(eq[0])
# 将x的解代入y的方程再求解
y_eq = eq[1].subs(x(t), x_sol.rhs)
y_sol = sp.dsolve(y_eq)

print(x_sol, y_sol)

运行后x(t)直接输出为Eq(x(t), C1*exp(-a*t))，形式和单方程求解完全一致，y(t)化简后和Mathematica结果完全匹配。

• 常数替换法：对Sympy返回的联立解做常数重定义，消去x(t)里的冗余参数

sol = sp.dsolve(eq)
# 提取解中的原有任意常数
old_consts = sorted([s for s in sol[0].rhs.free_symbols.union(sol[1].rhs.free_symbols) if str(s).startswith("C")], key=str)
C1, C2 = sp.symbols("C1 C2")
# 定义新旧常数的替换关系，让x(t)的系数归一
replace_map = {
    old_consts[0]: -a*C1/(a-b),
    old_consts[1]: C2 - a*C1/(a-b)
}
simplified_sol = [eq.subs(replace_map).simplify() for eq in sol]
print(simplified_sol)

提示：微分方程的通解中，任意常数的线性组合、缩放都不改变解的正确性，不同求解工具、不同求解顺序返回的解形式可能不同，只要代入原方程等式成立就是有效解。

内容的提问来源于stack exchange，提问作者BBQuercus