Sympy - ограничения в вычислительной способности функции dsolve?

Question

Я пытаюсь решить систему дифференциальных уравнений четвертого порядка, используя Sympy.Уравнения, которые я использовал, показаны на рисунке и написаны в следующем коде:

latex_equations:

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 больше не сможет ее решить?Любая помощь или понимание этой проблемы будет принята с благодарностью.

Спасибо!

Answer 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 в первом уравнении нетак же, как С2 во втором.Система ODE должна иметь 12 различных констант интеграции.На самом деле, хотя это система DAE, тогда должно быть 8 независимых констант, поскольку между этими тремя зависимыми переменными существует алгебраическая связь.