Search code examples
pythonsympydsolve

Sympy - limitations in the computation ability of dsolve function?

I am trying to solve a system of fourth-order differential equations using Sympy. The equations I have used are shown in the image, and written in the code below:

latex_equations :

enter image description here

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)

And I received the following error:

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'

I tried to solve for a simpler system of equations using dsolve, and it solved fine:

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)

The format I used for these two cases is the same, yet one solves and one fails. I know that the first system of equations has a solution, because I have solved for it in Maple.

Have I made an error in my code, or is Sympy dsolve simply not able to solve such a complex system of equations? Is there a limitation to how complex the system of equations may be until dsolve can no longer solve it? Any help or insight into this problem would greatly be appreciated.

Thank you!

Solution

  • SymPy's ODE systems code needs a lot of work and can't currently handle this kind of system.

    Firstly the third equation is an algebraic equation (meaning that this is really a DAE rather than ODE system) but that's not a huge problem: just differentiate it 4 times. Then SymPy should really be able to solve the resulting system but unfortunately can't because it's not implemented yet. You can do it yourself manually though:

    # 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))

    This gives the solutions:

                                                         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₃)

    Note though that the integration constants for each should be different so C1 in the first equation is not the same as C2 in the second. The ODE system should have 12 different integration constants. Actually though since this is a DAE system then there should be 8 independent constants since there is an algebraic relationship between the three dependent variables.

    EDIT: In SymPy 1.8 the dsolve function can now solve this provided the 3rd equation is differentiated:

    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₃))⎦