在 python (sympy) 中使用 dsolve_system 求解耦合微分方程组

如何解决在 python (sympy) 中使用 dsolve_system 求解耦合微分方程组

我想用 python (sympy) 求解一个由 4 个耦合的微分方程组成的系统：

eqs = [Eq(cP1(t).diff(t),k1*cE1(t)**3),Eq(cE1(t).diff(t),-k1 * cE1(t)**3 + k6 * cE3(t)**2),Eq(cE2(t).diff(t),-k8 * cE2(t)),Eq(cE3(t).diff(t),k8 * cE2(t) - k6 * cE3(t)**2)]

当我尝试用“dsolve_system”解决系统问题时：

solution = dsolve_system(eqs,ics={cP1(0): 0,cE1(0): cE1_0,cE2(0):cE2_0,cE3(0):cE3_0})

我得到答案：“NotImplementedError：
传递的 ODE 的系统无法通过 dsolve_system 求解。"

有谁知道，有什么问题吗？或者有没有更好的方法可以在 Sympy 中解决这个微分方程组？
非常感谢！

解决方法

显示完整代码通常很好且友好：

In [18]: cP1,cE1,cE2,cE3 = symbols('cP1,cE1:4',cls=Function)

In [19]: t,k1,k6,k8 = symbols('t,k8')

In [20]: eqs = [Eq(cP1(t).diff(t),k1*cE1(t)**3),Eq(cE1(t).diff(t),-k1 * cE1(t)**3 + k6 * cE3(t)**2),...:        Eq(cE2(t).diff(t),-k8 * cE2(t)),Eq(cE3(t).diff(t),k8 * cE2(t) - k6 * cE3(t)**2)]
    ...: 

In [21]: for eq in eqs:
    ...:     pprint(eq)
    ...: 
d                  3   
──(cP₁(t)) = k₁⋅cE₁ (t)
dt                     
d                    3            2   
──(cE₁(t)) = - k₁⋅cE₁ (t) + k₆⋅cE₃ (t)
dt                                    
d                      
──(cE₂(t)) = -k₈⋅cE₂(t)
dt                     
d                    2               
──(cE₃(t)) = - k₆⋅cE₃ (t) + k₈⋅cE₂(t)
dt                                   

In [22]: dsolve(eqs)
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-22-69ab769a7261> in <module>
----> 1 dsolve(eqs)

~/current/sympy/sympy/sympy/solvers/ode/ode.py in dsolve(eq,func,hint,simplify,ics,xi,eta,x0,n,**kwargs)
    609             "number of functions being equal to number of equations")
    610         if match['type_of_equation'] is None:
--> 611             raise NotImplementedError
    612         else:
    613             if match['is_linear'] == True:

NotImplementedError: 

这意味着 dsolve 还不能处理这种特定类型的系统。请注意，通常 ODE 的非线性系统不太可能有解析解（dsolve 用于寻找解析解，如果您想要数值解，请使用 scipy 的 odeint 之类的东西）。

随着非线性系统的发展，这个相对友好，所以有可能解决它。让我们看看...

首先，我们可以使用一个守恒量（所有 4 个变量的总和）来消除一个方程。实际上这并没有多大帮助，因为第一个方程已经与其他方程分离了：如果我们知道 cE1，我们就可以积分，但如果其他变量已知，守恒量就更容易给出它。

系统结构如下：

cE2   --->   cE3   --->   cE1  --->  cP1

这意味着它可以作为一系列 ODE 求解，其中我们求解 cE2 的第三个方程，然后求解 cE3 的第四个方程，然后将其用于 cE1，依此类推。因此，我们可以将其简化为涉及一系列主要是非线性的单个 ODE 的问题。这也不太可能有分析解决方案，但让我们试一试：

In [24]: dsolve(eqs[2],cE2(t))
Out[24]: 
             -k₈⋅t
cE₂(t) = C₁⋅ℯ     

In [25]: cE2sol = dsolve(eqs[2],cE2(t)).rhs

In [26]: eqs[3].subs(cE2(t),cE2sol)
Out[26]: 
d                   -k₈⋅t         2   
──(cE₃(t)) = C₁⋅k₈⋅ℯ      - k₆⋅cE₃ (t)
dt   

在这一点上，原则上我们可以在这里求解 cE3，但 sympy 没有任何方法来求解这个特定的非线性 ODE，所以 dsolve 给出了一个级数解（我认为这不是您想要），而唯一可以处理此问题的其他求解器是 lie_group，但实际上失败了。

由于我们无法得到 cE3 的解的表达式，我们也无法继续求解 cE1 和 cP1。失败的 ODE 有一个 Riccati 方程，但这种特殊类型的 Ricatti 方程尚未在 dsolve 中实现。看起来 Wolfram Alpha 给出了贝塞尔函数的答案： https://www.wolframalpha.com/input/?i=solve+dx%2Fdt+%3D+e%5E-t+-+x%5E2

从这一点来看，我猜我们不太可能解出下一个方程。那时 Wolfram Alpha 也放弃了并说“分类：非线性微分方程组”： https://www.wolframalpha.com/input/?i=solve+dx%2Fdt+%3D+e%5E-t+-+x%5E2%2C+dy%2Fdt+%3D+-y%5E3+%2B+x%5E2

我怀疑您的系统没有解析解，但您可以尝试数值解或其他更定性的分析。