0

我正在尝试使用 Sympy 求解四阶微分方程组。我使用的方程式如图所示，并写在下面的代码中：

乳胶方程：

在此处输入图像描述

from sympy import *
x = symbols('x')
EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
u1,u2,u3 = symbols('u1 u2 u3', cls=Function)

eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))

dsolve(eq)

我收到以下错误：

AttributeError                            Traceback (most recent call last)
<ipython-input-1-63c42d2751be> in <module>
      6 eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
      7 
----> 8 dsolve(eq)

~\Anaconda3\lib\site-packages\sympy\solvers\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
    583     """
    584     if iterable(eq):
--> 585         match = classify_sysode(eq, func)
    586         eq = match['eq']
    587         order = match['order']

~\Anaconda3\lib\site-packages\sympy\solvers\ode.py in classify_sysode(eq, funcs, **kwargs)
   1528         if isinstance(func, list):
   1529             for func_elem in func:
-> 1530                 if len(func_elem.args) != 1:
   1531                     raise ValueError("dsolve() and classify_sysode() work with "
   1532                     "functions of one variable only, not %s" % func)

AttributeError: 'list' object has no attribute 'args'

我尝试使用 dsolve 求解一个更简单的方程组，结果很好：

from sympy import *
x = symbols('x')
EI1,EI2 = symbols('EI1 EI2')
u1,u2 = symbols('u1 u2', cls=Function)

eq = (Eq(EI1*diff(u1(x),x), 12*x*u1(x) + 8*u2(x)), Eq(EI2*diff(u2(x),x), 21*u1(x) + 7*x*u2(x)))

dsolve(eq)

我在这两种情况下使用的格式是相同的，但一种解决了，一种失败了。我知道第一个方程组有一个解决方案，因为我已经在 Maple 中解决了它。

我的代码是否出错，或者 Sympy dsolve 根本无法解决如此复杂的方程组？在 dsolve 无法解决之前，方程组的复杂程度是否存在限制？对此问题的任何帮助或见解将不胜感激。

谢谢！

4

1 回答 1

1

SymPy 的 ODE 系统代码需要大量工作，目前无法处理这种系统。

首先，第三个方程是一个代数方程（这意味着这实际上是一个 DAE 而不是 ODE 系统）但这不是一个大问题：只需将其微分 4 次。那么 SymPy 应该真的能够解决生成的系统，但不幸的是不能，因为它还没有实现。您可以自己手动完成：

# Your code first:
from sympy import *
x = symbols('x')
EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
u1,u2,u3 = symbols('u1 u2 u3', cls=Function)

eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3*diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))

# Differentiate eq3 4 times:
eq1, eq2, eq3 = eq
eq3 = Eq(eq3.lhs.diff(x, 4).doit(), eq3.rhs.diff(x, 4).doit())
eq = eq1, eq2, eq3

# Isolate equations for u1(x).diff(x, 4) etc
derivs = [u(x).diff(x, 4) for u in (u1, u2, u3)]
(sol,) = solve(eq, derivs, dict=True)
eq = [Eq(du, sol[du]) for du in derivs]

# Each can be solved separately:
for e in eq:
    pprint(dsolve(e))

这给出了解决方案：

                                                     4                                                                  
                        2       3                   x ⋅(EI₂⋅Mecc - EI₂⋅Qh⋅a₂ + 2⋅EI₃⋅Mecc - 2⋅EI₃⋅Qh⋅a₃)                
u₁(x) = C₁ + C₂⋅x + C₃⋅x  + C₄⋅x  + ────────────────────────────────────────────────────────────────────────────────────
                                    24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)
                                                      4                                                                 
                        2       3                    x ⋅(-EI₁⋅Mecc + EI₁⋅Qh⋅a₁ + EI₃⋅Mecc - EI₃⋅Qh⋅a₃)                  
u₂(x) = C₁ + C₂⋅x + C₃⋅x  + C₄⋅x  + ────────────────────────────────────────────────────────────────────────────────────
                                    24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)
                                                    4                                                                   
                        2       3                  x ⋅(-2⋅EI₁⋅Mecc + 2⋅EI₁⋅Qh⋅a₁ - EI₂⋅Mecc + EI₂⋅Qh⋅a₂)                
u₃(x) = C₁ + C₂⋅x + C₃⋅x  + C₄⋅x  + ────────────────────────────────────────────────────────────────────────────────────
                                    24⋅(EI₁⋅EI₂⋅a₁ - EI₁⋅EI₂⋅a₂ + 2⋅EI₁⋅EI₃⋅a₁ - 2⋅EI₁⋅EI₃⋅a₃ + EI₂⋅EI₃⋅a₂ - EI₂⋅EI₃⋅a₃)

请注意，每个方程的积分常数应该不同，因此第一个方程中的 C1 与第二个方程中的 C2 不同。ODE 系统应该有 12 个不同的积分常数。实际上，由于这是一个 DAE 系统，因此应该有 8 个独立常数，因为三个因变量之间存在代数关系。

编辑：在 SymPy 1.8 中，dsolve函数现在可以解决这个问题，前提是第三个方程是微分的：

In [67]: # Your code first:
    ...: from sympy import *
    ...: x = symbols('x')
    ...: EI1,EI2,EI3,a1,a2,a3,Qh,Mecc = symbols('EI1 EI2 EI3 a1 a2 a3 Qh Mecc')
    ...: u1,u2,u3 = symbols('u1 u2 u3', cls=Function)
    ...: 
    ...: eq = (Eq(EI1*diff(u1(x),x,x,x,x)+EI2*diff(u2(x),x,x,x,x)+EI3*diff(u3(x),x,x,x,x), Qh),Eq(a1*EI1*diff(u1(x),x,x,x,x)+a2*EI2*diff(u2(x),x,x,x,x)+a3*EI3
    ...: *diff(u3(x),x,x,x,x),Mecc),Eq((u1(x)+u3(x))/2,u2(x)))
    ...: 
    ...: # Differentiate eq3 4 times:
    ...: eq1, eq2, eq3 = eq
    ...: eq3 = Eq(eq3.lhs.diff(x, 4).doit(), eq3.rhs.diff(x, 4).doit())
    ...: eq = eq1, eq2, eq3

In [68]: dsolve(eq)
Out[68]: 
⎡                        2       3            4                                                                              2       3             4        
⎢                    C₃⋅x    C₄⋅x            x ⋅(EI₂⋅(Mecc - Qh⋅a₂) + 2⋅EI₃⋅(Mecc - Qh⋅a₃))                              C₇⋅x    C₈⋅x             x ⋅(EI₁⋅(M
⎢u₁(x) = C₁ + C₂⋅x + ───── + ───── + ──────────────────────────────────────────────────────────────, u₂(x) = C₅ + C₆⋅x + ───── + ───── - ───────────────────
⎣                      2       6     24⋅(EI₁⋅(EI₂⋅(a₁ - a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃))                        2       6     24⋅(EI₁⋅(EI₂⋅(a₁ - 

                                                                  2        3                 4                                                    ⎤
ecc - Qh⋅a₁) - EI₃⋅(Mecc - Qh⋅a₃))                           C₁₁⋅x    C₁₂⋅x                 x ⋅(2⋅EI₁⋅(Mecc - Qh⋅a₁) + EI₂⋅(Mecc - Qh⋅a₂))        ⎥
───────────────────────────────────────────, u₃(x) = C₁₀⋅x + ────── + ────── + C₉ - ──────────────────────────────────────────────────────────────⎥
a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃))                    2        6           24⋅(EI₁⋅(EI₂⋅(a₁ - a₂) + 2⋅EI₃⋅(a₁ - a₃)) + EI₂⋅EI₃⋅(a₂ - a₃))⎦
于 2019-05-22T22:09:24.263 回答