0

我想用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 中解决这个微分方程系统？
非常感谢！

4

1 回答 1

0

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

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

In [19]: t, k1, k6, k8 = symbols('t, k1, k6, 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

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

于 2020-12-21T15:12:21.190 回答