New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

dsolve system gives complicated solution for diagonal system #15474

Open
oscarbenjamin opened this Issue Nov 9, 2018 · 3 comments

Comments

Projects
None yet
2 participants
@oscarbenjamin
Copy link
Contributor

oscarbenjamin commented Nov 9, 2018

This is a very simple system of ODEs:

In [1]: x, y = symbols('x y', cls=Function)

In [2]: a, b, t = symbols('a b t')

In [3]: eq = [Eq(x(t).diff(t), a*x(t)), Eq(y(t).diff(t), b*y(t))]

In [5]: eq
Out[5]: 
⎡d                  d                ⎤
⎢──(x(t)) = a⋅x(t), ──(y(t)) = b⋅y(t)⎥
⎣dt                 dt               ⎦

The solution from sympy explodes (see below). The correct solution is

In [9]: [Eq(x(t), C1*exp(a*t)), Eq(y(t), C2*exp(b*t))]
Out[9]: 
⎡           a⋅t             b⋅t⎤
⎣x(t) = C₁⋅ℯ   , y(t) = C₂⋅ℯ   ⎦

Calling dsolve here calls through to _linear_2eq_order1_type1. Having looked at this function I don't understand why it is written the way that it is with so many piecewise parts to the solution (I do understand why they are technically correct but I don't think it's reasonable to try to catch all cases with piecewise like this). The solution to the system mentioned here is easily found using _linear_neq_order1_type1 which should in general make _linear_2eq_order1_type1 redundant. I think that the best fix here is simply to remove _linear_2eq_order1_type1.

I just tested this diff and all tests pass except for test_dsolve_linsystem_symbol_piecewise which specifically tests for the kind of complicated piecewise expression seen in the undesired output below:

diff --git a/sympy/solvers/ode.py b/sympy/solvers/ode.py
index 8fc6f514a..cb0657b97 100644
--- a/sympy/solvers/ode.py
+++ b/sympy/solvers/ode.py
@@ -6745,7 +6745,7 @@ def sysode_linear_2eq_order1(match_):
                                   " (and constant inhomogeneity)")
 
     if match_['type_of_equation'] == 'type1':
-        sol = _linear_2eq_order1_type1(x, y, t, r, eq)
+        sol = _linear_neq_order1_type1(match_)

Here is the full output from current master for this system of ODEs:

In [4]: dsolve(eq)
Out[4]: 
⎡                                                                                                                                            
⎢                                                                                                                                            
⎢       ⎧                                                                      ⎛           __________⎞                                       
⎢       ⎪                                                                      ⎜          ╱        2 ⎟                                       
⎢       ⎪           ⎛⎧           __________                     __________⎞    ⎜a   b   ╲╱  (a - b)  ⎟      ⎛⎧           __________          
⎢       ⎪           ⎜⎪          ╱        2                     ╱        2 ⎟  t⋅⎜─ + ─ + ─────────────⎟      ⎜⎪          ╱        2           
⎢       ⎪           ⎜⎪a   b   ╲╱  (a - b)            a   b   ╲╱  (a - b)  ⎟    ⎝2   2         2      ⎠      ⎜⎪a   b   ╲╱  (a - b)            
⎢       ⎪        C₁⋅⎜⎨─ - ─ + ─────────────  for a = ─ + ─ + ─────────────⎟⋅ℯ                          + C₂⋅⎜⎨─ - ─ - ─────────────  for a = 
⎢       ⎪           ⎜⎪2   2         2                2   2         2      ⎟                                 ⎜⎪2   2         2                
⎢       ⎪           ⎜⎪                                                    ⎟                                 ⎜⎪                               
⎢       ⎪           ⎝⎩          0                      otherwise          ⎠                                 ⎝⎩          0                    
⎢       ⎪                                                                                                                                    
⎢       ⎪                                                                                      ⎛           __________⎞                       
⎢       ⎪                                                                                      ⎜          ╱        2 ⎟                       
⎢       ⎪                                                                                      ⎜a   b   ╲╱  (a - b)  ⎟                       
⎢x(t) = ⎨                                                                                    t⋅⎜─ + ─ + ─────────────⎟                       
⎢       ⎪                                                                                      ⎝2   2         2      ⎠                       
⎢       ⎪                                                                                C₁⋅ℯ                                                
⎢       ⎪                                                                                                                                    
⎢       ⎪                                                                                                                                    
⎢       ⎪                                                                                                                                    
⎢       ⎪⎛   ⎛⎧           __________                     __________⎞      ⎛  ⎛⎧           __________                     __________⎞   ⎛⎧    
⎢       ⎪⎜   ⎜⎪          ╱        2                     ╱        2 ⎟      ⎜  ⎜⎪          ╱        2                     ╱        2 ⎟   ⎜⎪    
⎢       ⎪⎜   ⎜⎪a   b   ╲╱  (a - b)            a   b   ╲╱  (a - b)  ⎟      ⎜  ⎜⎪a   b   ╲╱  (a - b)            a   b   ╲╱  (a - b)  ⎟   ⎜⎪    
⎢       ⎪⎜C₁⋅⎜⎨─ - ─ + ─────────────  for a = ─ + ─ + ─────────────⎟ + C₂⋅⎜t⋅⎜⎨─ - ─ + ─────────────  for a = ─ + ─ + ─────────────⎟ + ⎜⎨1  f
⎢       ⎪⎜   ⎜⎪2   2         2                2   2         2      ⎟      ⎜  ⎜⎪2   2         2                2   2         2      ⎟   ⎜⎪    
⎢       ⎪⎜   ⎜⎪                                                    ⎟      ⎜  ⎜⎪                                                    ⎟   ⎜⎪    
⎢       ⎪⎝   ⎝⎩          0                      otherwise          ⎠      ⎝  ⎝⎩          0                      otherwise          ⎠   ⎝⎩0   
⎢       ⎩                                                                                                                                    
⎢                                                                                                                                            
⎣                                                                                                                                            

                                                                                            ⎧   ⎛⎧                                           
                                                                                            ⎪   ⎜⎪                                           
                          ⎛           __________⎞                                           ⎪   ⎜⎪                                 a   b   ╲╱
                          ⎜          ╱        2 ⎟                                           ⎪   ⎜⎪           0             for a = ─ + ─ + ──
           __________⎞    ⎜a   b   ╲╱  (a - b)  ⎟                                           ⎪   ⎜⎪                                 2   2     
          ╱        2 ⎟  t⋅⎜─ + ─ - ─────────────⎟                                           ⎪C₁⋅⎜⎨                                           
a   b   ╲╱  (a - b)  ⎟    ⎝2   2         2      ⎠                2            2             ⎪   ⎜⎪             __________                    
─ + ─ - ─────────────⎟⋅ℯ                                    for a  - 2⋅a⋅b + b  ≠ 0         ⎪   ⎜⎪            ╱        2                     
2   2         2      ⎟                                                                      ⎪   ⎜⎪  a   b   ╲╱  (a - b)                      
                     ⎟                                                                      ⎪   ⎜⎪- ─ + ─ + ─────────────            otherwis
  otherwise          ⎠                                                                      ⎪   ⎝⎩  2   2         2                          
                                                                                            ⎪                                                
                                                                                            ⎪                                                
                                                                                            ⎪                                                
                                                                                            ⎪                                                
                                                                                   , y(t) = ⎨                                                
                                                                                            ⎪                                                
                                                                   for a = b                ⎪                                                
                                                                                            ⎪                                                
                                   ⎛           __________⎞                                  ⎪           ⎛   ⎛⎧                               
                                   ⎜          ╱        2 ⎟                                  ⎪           ⎜   ⎜⎪                               
                  __________⎞⎞⎞    ⎜a   b   ╲╱  (a - b)  ⎟                                  ⎪           ⎜   ⎜⎪                               
                 ╱        2 ⎟⎟⎟  t⋅⎜─ + ─ + ─────────────⎟                                  ⎪           ⎜   ⎜⎪           0             for a 
       a   b   ╲╱  (a - b)  ⎟⎟⎟    ⎝2   2         2      ⎠                                  ⎪           ⎜   ⎜⎪                               
or a = ─ + ─ + ─────────────⎟⎟⎟⋅ℯ                                  otherwise                ⎪           ⎜C₁⋅⎜⎨                               
       2   2         2      ⎟⎟⎟                                                             ⎪           ⎜   ⎜⎪             __________        
                            ⎟⎟⎟                                                             ⎪           ⎜   ⎜⎪            ╱        2         
         otherwise          ⎠⎠⎠                                                             ⎪           ⎜   ⎜⎪  a   b   ╲╱  (a - b)          
                                                                                            ⎪           ⎜   ⎜⎪- ─ + ─ + ─────────────        
                                                                                            ⎪           ⎝   ⎝⎩  2   2         2              
                                                                                            ⎩                                                

 __________⎞    ⎛           __________⎞      ⎛⎧                                            __________⎞    ⎛           __________⎞            
╱        2 ⎟    ⎜          ╱        2 ⎟      ⎜⎪                                           ╱        2 ⎟    ⎜          ╱        2 ⎟            
  (a - b)  ⎟    ⎜a   b   ╲╱  (a - b)  ⎟      ⎜⎪                                 a   b   ╲╱  (a - b)  ⎟    ⎜a   b   ╲╱  (a - b)  ⎟            
───────────⎟  t⋅⎜─ + ─ + ─────────────⎟      ⎜⎪           0             for a = ─ + ─ - ─────────────⎟  t⋅⎜─ + ─ - ─────────────⎟            
    2      ⎟    ⎝2   2         2      ⎠      ⎜⎪                                 2   2         2      ⎟    ⎝2   2         2      ⎠       2    
           ⎟⋅ℯ                          + C₂⋅⎜⎨                                                      ⎟⋅ℯ                           for a  - 2
           ⎟                                 ⎜⎪             __________                               ⎟                                       
           ⎟                                 ⎜⎪            ╱        2                                ⎟                                       
           ⎟                                 ⎜⎪  a   b   ╲╱  (a - b)                                 ⎟                                       
e          ⎟                                 ⎜⎪- ─ + ─ - ─────────────            otherwise          ⎟                                       
           ⎠                                 ⎝⎩  2   2         2                                     ⎠                                       
                                                                                                                                             
                                ⎛           __________⎞                                                                                      
                                ⎜          ╱        2 ⎟                                                                                      
                                ⎜a   b   ╲╱  (a - b)  ⎟                                                                                      
                              t⋅⎜─ + ─ + ─────────────⎟                                                                                      
                                ⎝2   2         2      ⎠                                                                                      
                          C₂⋅ℯ                                                                                                            for
                                                                                                                                             
             __________⎞        ⎛⎧                                            __________⎞⎞    ⎛           __________⎞                        
            ╱        2 ⎟        ⎜⎪                                           ╱        2 ⎟⎟    ⎜          ╱        2 ⎟                        
  a   b   ╲╱  (a - b)  ⎟        ⎜⎪                                 a   b   ╲╱  (a - b)  ⎟⎟    ⎜a   b   ╲╱  (a - b)  ⎟                        
= ─ + ─ + ─────────────⎟        ⎜⎪           0             for a = ─ + ─ + ─────────────⎟⎟  t⋅⎜─ + ─ + ─────────────⎟                        
  2   2         2      ⎟        ⎜⎪                                 2   2         2      ⎟⎟    ⎝2   2         2      ⎠                        
                       ⎟ + C₂⋅t⋅⎜⎨                                                      ⎟⎟⋅ℯ                                              oth
                       ⎟        ⎜⎪             __________                               ⎟⎟                                                   
                       ⎟        ⎜⎪            ╱        2                                ⎟⎟                                                   
                       ⎟        ⎜⎪  a   b   ╲╱  (a - b)                                 ⎟⎟                                                   
    otherwise          ⎟        ⎜⎪- ─ + ─ + ─────────────            otherwise          ⎟⎟                                                   
                       ⎠        ⎝⎩  2   2         2                                     ⎠⎠                                                   
                                                                                                                                             

             ⎤
             ⎥
             ⎥
             ⎥
        2    ⎥
⋅a⋅b + b  ≠ 0⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
 a = b       ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
erwise       ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎥
             ⎦

@oscarbenjamin oscarbenjamin changed the title dsolve system gives complex solution for diagonal system dsolve system gives complicated solution for diagonal system Nov 10, 2018

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 10, 2018

CC @cbm755 who may have opinions no the systems solver.

skirpichev added a commit to skirpichev/diofant that referenced this issue Nov 10, 2018

Add regression test
Closes sympy/sympy#15474

N.b., in v0.7.6 this raises a PolynomialError; SymPy got a
complicated answer instead in v1.2.
@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

oscarbenjamin commented Nov 10, 2018

I have another example:

In [1]: x, y = symbols('x y', cls=Function)

In [2]: a = Symbol('a')

In [3]: eq = [Eq(x(t).diff(t), -a*y(t)), Eq(y(t).diff(t), a*x(t))]

In [4]: eq
Out[4]: 
⎡d                   d                ⎤
⎢──(x(t)) = -a⋅y(t), ──(y(t)) = a⋅x(t)⎥
⎣dt                  dt               ⎦

The correct solution is

In [8]: [Eq(x(t), C1*cos(a*t) + C2*sin(a*t)), Eq(y(t), C1*sin(a*t) - C2*cos(a*t))]
Out[8]: [x(t) = C₁⋅cos(a⋅t) + C₂⋅sin(a⋅t), y(t) = C₁⋅sin(a⋅t) - C₂⋅cos(a⋅t)]

Full output from sympy:

In [5]: dsolve(eq)
Out[5]: 
⎡       ⎧                                            _____                                     _____                                         
⎢       ⎪               ⎛⎧   _____           ⎞      ╱   2       ⎛⎧    _____           ⎞       ╱   2                                          
⎢       ⎪               ⎜⎪  ╱   2            ⎟  t⋅╲╱  -a        ⎜⎪   ╱   2-t⋅╲╱  -a                        2             ⎧    
⎢       ⎪            C₁⋅⎜⎨╲╱  -a    for a = 0⎟⋅ℯ           + C₂⋅⎜⎨-╲╱  -a    for a = 0⎟⋅ℯ                         for -4⋅a  ≠ 0         ⎪    
⎢       ⎪               ⎜⎪                   ⎟                  ⎜⎪                    ⎟                                                 ⎪    
⎢       ⎪               ⎝⎩   -a     otherwise⎠                  ⎝⎩   -a      otherwise⎠                                                 ⎪    
⎢       ⎪                                                                                                                               ⎪    
⎢       ⎪                                                     _____                                                                     ⎪    
⎢       ⎪                                                    ╱   2                                                                      ⎪    
⎢       ⎪                                                t⋅╲╱  -a                                                                       ⎪    
⎢       ⎪                                            C₁⋅ℯ                                                           for a = 0           ⎪    
⎢x(t) = ⎨                                                                                                                      , y(t) = ⎨    
⎢       ⎪⎛                               ⎛                           ⎛⎧   1         for a = 0    ⎞⎞⎞                                    ⎪    
⎢       ⎪⎜                               ⎜                           ⎜⎪                          ⎟⎟⎟                                    ⎪    
⎢       ⎪⎜                               ⎜                           ⎜⎪                 _____    ⎟⎟⎟       _____                        ⎪⎛   
⎢       ⎪⎜   ⎛⎧   _____           ⎞      ⎜  ⎛⎧   _____           ⎞   ⎜⎪                ╱   2     ⎟⎟⎟      ╱   2                         ⎪⎜   
⎢       ⎪⎜   ⎜⎪  ╱   2            ⎟      ⎜  ⎜⎪  ╱   2            ⎟   ⎜⎪   0      for ╲╱  -a   = 0⎟⎟⎟  t⋅╲╱  -a                          ⎪⎜   
⎢       ⎪⎜C₁⋅⎜⎨╲╱  -a    for a = 0+ C₂⋅⎜t⋅⎜⎨╲╱  -a    for a = 0+ ⎜⎨                          ⎟⎟⎟⋅ℯ              otherwise           ⎪⎜C₁⋅
⎢       ⎪⎜   ⎜⎪                   ⎟      ⎜  ⎜⎪                   ⎟   ⎜⎪   a                      ⎟⎟⎟                                    ⎪⎜   
⎢       ⎪⎜   ⎝⎩   -a     otherwise⎠      ⎜  ⎝⎩   -a     otherwise⎠   ⎜⎪────────     otherwise    ⎟⎟⎟                                    ⎪⎜   
⎢       ⎪⎜                               ⎜                           ⎜⎪   _____                  ⎟⎟⎟                                    ⎩⎝   
⎢       ⎪⎜                               ⎜                           ⎜⎪  ╱   2                   ⎟⎟⎟                                         
⎣       ⎩⎝                               ⎝                           ⎝⎩╲╱  -a                    ⎠⎠⎠                                         

                                                                                                           ⎤
                                                                                                           ⎥
                                    _____                                     _____                        ⎥
       ⎛⎧   a      for a = 0⎞      ╱   2       ⎛⎧    a      for a = 0⎞       ╱   2                         ⎥
       ⎜⎪                   ⎟  t⋅╲╱  -a        ⎜⎪                    ⎟  -t⋅╲╱  -a                     2    ⎥
    C₁⋅⎜⎨   _____           ⎟⋅ℯ           + C₂⋅⎜⎨    _____           ⎟⋅ℯ                      for -4⋅a  ≠ 0⎥
       ⎜⎪  ╱   2            ⎟                  ⎜⎪   ╱   2            ⎟                                     ⎥
       ⎝⎩╲╱  -a    otherwise⎠                  ⎝⎩-╲╱  -a    otherwise⎠                                     ⎥
                                                                                                           ⎥
                                              _____                                                        ⎥
                                             ╱   2                                                         ⎥
                                         t⋅╲╱  -a                                                          ⎥
                                     C₂⋅ℯ                                                       for a = 0  ⎥
                                                                                                           ⎥
                            ⎛                           ⎛⎧0     for a = 0    ⎞⎞⎞       _____               ⎥
⎛⎧   a      for a = 0⎞      ⎜  ⎛⎧   a      for a = 0⎞   ⎜⎪                   ⎟⎟⎟      ╱   2                ⎥
⎜⎪                   ⎟      ⎜  ⎜⎪                   ⎟   ⎜⎪          _____    ⎟⎟⎟  t⋅╲╱  -a                 ⎥
⎜⎨   _____           ⎟ + C₂⋅⎜t⋅⎜⎨   _____           ⎟ + ⎜⎨         ╱   2     ⎟⎟⎟⋅ℯ              otherwise  ⎥
⎜⎪  ╱   2            ⎟      ⎜  ⎜⎪  ╱   2            ⎟   ⎜⎪1  for ╲╱  -a   = 0⎟⎟⎟                           ⎥
⎝⎩╲╱  -a    otherwise⎠      ⎜  ⎝⎩╲╱  -a    otherwise⎠   ⎜⎪                   ⎟⎟⎟                           ⎥
                            ⎝                           ⎝⎩0     otherwise    ⎠⎠⎠                           ⎥
                                                                                                           ⎥
                                                                                                           ⎦

And the same with a declared as real:

In [6]: a = Symbol('a', real=True)

In [7]: eq = [Eq(x(t).diff(t), -a*y(t)), Eq(y(t).diff(t), a*x(t))]

In [8]: dsolve(eq)
Out[8]: 
⎡       ⎧                                                                                                  2                                 
⎢       ⎪                           -C₁⋅a⋅cos(t⋅│a│) + C₂⋅a⋅sin(t⋅│a│)                             for -4⋅a  < 0         ⎧                   
⎢       ⎪                                                                                                                ⎪                   
⎢       ⎪                                          ⅈ⋅t⋅│a│                                                               ⎪                   
⎢       ⎪                                      C₁⋅ℯ                                                  for a = 0           ⎪                   
⎢       ⎪                                                                                                                ⎪                   
⎢x(t) = ⎨⎛                            ⎛                        ⎛⎧  1      for a = 0  ⎞⎞⎞                        , y(t) = ⎨                   
⎢       ⎪⎜                            ⎜                        ⎜⎪                    ⎟⎟⎟                                 ⎪⎛                  
⎢       ⎪⎜   ⎛⎧ⅈ⋅│a│  for a = 0⎞      ⎜  ⎛⎧ⅈ⋅│a│  for a = 0⎞   ⎜⎪  0    for ⅈ⋅│a│ = 0⎟⎟⎟  ⅈ⋅t⋅│a│          2             ⎪⎜   ⎛⎧  a    for a 
⎢       ⎪⎜C₁⋅⎜⎨                ⎟ + C₂⋅⎜t⋅⎜⎨                ⎟ + ⎜⎨                    ⎟⎟⎟⋅ℯ         for -4⋅a  = 0         ⎪⎜C₁⋅⎜⎨             
⎢       ⎪⎜   ⎝⎩ -a    otherwise⎠      ⎜  ⎝⎩ -a    otherwise⎠   ⎜⎪-ⅈ⋅a                ⎟⎟⎟                                 ⎪⎜   ⎝⎩ⅈ⋅│a│  otherw
⎢       ⎪⎜                            ⎜                        ⎜⎪─────    otherwise  ⎟⎟⎟                                 ⎩⎝                  
⎣       ⎩⎝                            ⎝                        ⎝⎩ │a│                ⎠⎠⎠                                                     

                                                                                ⎤
                                                                           2-C₁⋅sin(t⋅│a│)⋅│a│ - C₂⋅cos(t⋅│a│)⋅│a│                         for -4⋅a  < 0⎥
                                                                                ⎥
                     ⅈ⋅t⋅│a│                                                    ⎥
                 C₂⋅ℯ                                                for a = 0  ⎥
                                                                                ⎥
          ⎛                        ⎛⎧0    for a = 0  ⎞⎞⎞                        ⎥
= 0⎞      ⎜  ⎛⎧  a    for a = 0⎞   ⎜⎪                ⎟⎟⎟  ⅈ⋅t⋅│a│          2    ⎥
   ⎟ + C₂⋅⎜t⋅⎜⎨                ⎟ + ⎜⎨1  for ⅈ⋅│a│ = 0⎟⎟⎟⋅ℯ         for -4⋅a  = 0⎥
ise⎠      ⎜  ⎝⎩ⅈ⋅│a│  otherwise⎠   ⎜⎪                ⎟⎟⎟                        ⎥
          ⎝                        ⎝⎩0    otherwise  ⎠⎠⎠                        ⎥
                                                                                ⎦

skirpichev added a commit to skirpichev/diofant that referenced this issue Nov 10, 2018

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

oscarbenjamin commented Nov 11, 2018

In fact there are many examples, I won't list them but I am adding them as tests to #15449 .

The piecewise code in _linear_2eq_order1_type1 is apparently added to handle the symbolic case but as far as I can tell it fails for almost all cases involving symbols.

You can see roughly what the output from _linear_neq_order1_type1 would be in each case by appending an additional trivial equation for z(t) e.g. for the first example in this issue:

In [1]: x, y, z = symbols('x y z', cls=Function)

In [2]: a, b, t = symbols('a b t')

In [3]: eq = [Eq(x(t).diff(t), a*x(t)), Eq(y(t).diff(t), b*y(t)), Eq(z(t).diff(t), 0)]

In [4]: dsolve(eq)
Out[4]: 
⎡           a⋅t             b⋅t           ⎤
⎣x(t) = C₁⋅ℯ   , y(t) = C₂⋅ℯ   , z(t) = C₃⎦

skirpichev added a commit to skirpichev/diofant that referenced this issue Nov 13, 2018

Add regression tests
Closes sympy/sympy#15474

N.b., in v0.7.6 this raises a PolynomialError; SymPy got a
complicated answer instead in v1.2.

skirpichev added a commit to skirpichev/diofant that referenced this issue Nov 13, 2018

Add regression tests
Closes sympy/sympy#15474

N.b., in v0.7.6 this raises a PolynomialError; SymPy got a
complicated answer instead in v1.2.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment