Tracking changes and tasks in systems.py

[mijo2]

This issue is created in order to keep track of the tasks that are needed to be done to extend the ODE module of SymPy. The solvers and methods to be added and that are already added are tracked below:

n-equations linear first order solvers:

• [x] constant coefficient homogeneous solver.

• [x] constant coefficient non-homogeneous solver.

• [x] non-constant coefficient homogeneous solver.

• [x] non-constant coefficient non-homogeneous solver.

Other functionalities:

• [x] Higher order to first order reduction.

• [x] Weakly and Strongly connected components.

• [ ] Special case non-linear solver where the independent variable can be eliminated.

• [ ] Undetermined coefficients solver.

But, before incorporating weakly and strongly connected components in solving the systems of ODEs, one of the main things that is needed is to add a main function named dsolve_system, which takes a system of ODEs as input and accordingly matches the system to one of the solvers or divides the system of ODEs if it is divisible using the concept of strongly and weakly connected components.

• [x] Adding the dsolve_system function along with the necessary helper functions.

Open PRs:

1. Simplification Strategies

Merged PRs:

1. Constant coefficient homogeneous solver
2. Non-constant coefficient homogeneous solver
3. Constant coefficient non-homogeneous solver
4. Non-constant coefficient non-homogeneous solver
5. API change for the newly added solvers
6. dsolve_system function
7. Tasks related to checksysodesol, constants_renumber and linear_ode_to_matrix
8. Strongly and weakly connected components to solve systems of ODEs
9. Higher-order to first-order reduction

More information about the tasks given in: https://github.com/sympy/sympy/wiki/ODE-Systems-roadmap

[mijo2]

I would say that the next step is to implement whichever part it is that can eliminate as much as possible of the old code. I'm not sure if that's the non-homogeneous case or the non-constant coefficient homogeneous case.

@oscarbenjamin I looked at the list of the solvers to be replaced when we add each of them. Adding non-constant coefficient homogeneous case would be wise since it replaces 10 solvers from the old code while, adding constant non-homogeneous case will result in removal of only one solver.

[oscarbenjamin]

Adding non-constant coefficient homogeneous case would be wise

I agree.

Another point to add to the todo list although it isn't in the roadmap is to make a systems version of the undetermined coefficients solver.

[sylee957]

How difficult is it to implement non-constant non-homogeneous (non-commutative) solver? I think that it can eventually be replaced with the most generalized solver. Are there any known answers?

[oscarbenjamin]

How difficult is it to implement non-constant non-homogeneous (non-commutative) solver?

In full generality this would imply a solver that is also capable of solving any single linear ODE including 2nd order examples like Bessel's equation, hypergeometric, Legendre etc. but also applying to arbitrary linear ODEs of any order. That does not seem easy to me!

There might be other general-ish cases that are solvable. The key question in terms of how solutions are presented in the roadmap is this:

• What you can say about d/dt exp(B(t)) based on knowing some properties of the matrix B?

So far all I have is based on B'(t) commuting with B(t).

[mijo2]

@oscarbenjamin I have updated the list since the second solver has been added. Now, should I start working on the constant coefficient non-homogeneous solver, since it will replace two old solvers along with the fact that it will lead to the closing of some of the issues related to solvers.dsolve.system

[oscarbenjamin]

Yes, that sounds good. That will be a significant extension of current capability so you should make sure to devise good tests for the new solver.

Ultimately we probably want to have an undetermined coefficients solver for many of those cases but since the integral solver is more general it's better to add that first.

[oscarbenjamin]

The constant coefficient non-homogeneous case is also needed by @namannimmo10 in the control theory GSOC project (#18436).

[oscarbenjamin]

The OP in this issue should be updated with references to PRs that are opened or merged. Then anyone looking at this can follow through to the in progress or completed PRs.

[moorepants]

What kinds of non-homogenous solutions to you all plan to support? There are a large set that would be useful for the controls module: various polynomials, sinusoidal, fourier series, etc.

[oscarbenjamin]

This is described in the ODE system roadmap: https://github.com/sympy/sympy/wiki/ODE-Systems-roadmap#constant-coefficient-non-homogeneous

The first step is to support these in integral form. Basically given

X' = AX + f(t)

with X and f vectors and A a matrix the solution is

X = exp(A t) int_x { exp(- A t) f(t) } + exp(A t) C

where C is a vector of integration constants.

This works for arbitrary non-homogeneous terms f(t) provided the integrals can be evaluated. Later it would probably be good to have a more efficient undetermined coefficients solver that would be applicable only to limited terms on the rhs (basically combinations of polynomial, exp, sin and cos).

[oscarbenjamin]

The other significant limitation that I should mention is needing to compute the matrix exponential which in turn requires being able to compute the eigenvalues.

[moorepants]

Sounds good. Yes, the eigenvalues will be an issue for symbolic forms and is of course limited by the size of A. I suspect it should return itself if it can't be solved and if it can, return the solution.

[oscarbenjamin]

There should probably be a way to have an unevaluated matrix exponential. We need a separate matrix exponential class rather than using matrices in the normal exp class.

[mijo2]

The constant coefficient non-homogeneous case is also needed by @namannimmo10 in the control theory GSOC project (#18436).

I will open a PR for the solver soon.

[mijo2]

@oscarbenjamin I will create new PR for the 4th linear solver soon.

[mijo2]

@oscarbenjamin As discussed, before moving on to the main function that will handle the solvers and the helper functions, I will first create a PR that will ensure that we can use the newly added solver without going through all the match functions and even an end-user can use the new solver by passing appropriate arguments.

[mijo2]

@oscarbenjamin I think before moving on to working with component division, it would be wise to make a PR that completes the points that you referenced in these comments:

1. Synchronising the outputs of linear_ode_to_matrix and canonical_odes
2. Adding checksysodesol call in checkodesol
3. Fixing constants_renumber function to handle list of ODEs

If there isn't an issue with creating a PR that addresses all the three comments, then I would create one soon.

[oscarbenjamin]

That seems reasonable to me.

[oscarbenjamin]

Another step to add is having a solver for Cauchy-Euler type systems like

t^2 A X'' + t B X' + C X = d(t)         (1)

Those can be solved with the substitution tau = ln(t) for the independent variable. Then if X(t) = Y(tau) = Y(ln(t)) we have

X'(t) = d/dt X(t) = d/dt Y(tau) = Y'(tau) dtau/dt = Y'(tau) / t

That means that t X' becomes Y' in the substitution. Likewise for X'' we have

X'' = d/dt (X') = d/dt (Y' / t) = Y'' / t^2 - Y' / t^2 = (Y'' - Y') / t^2

Now t^2 X'' becomesY'' - Y'. The same approach works for any order (not just 2nd order) and any number of equations and will turn the system into a constant coefficient system. It can then be solved astype1ortype2` depending on whether it is homogeneous or not.

Also 1st order systems of the form

X' = f(t) A X + b(t)            (2)

where f(t) is a scalar function and A is a constant matrix can be solved more simply using a subsitution tau = F(t) with X(t) = Y(tau) = Y(F(t)). We get

X' = d/dt X(t) = d/dt Y(tau) = Y'(tau) dtau/dt = Y'(tau) F'(t) = Y' F'(t)

so the transformed system is

Y' F'(t) = f(t) A Y + b(t)

If we choose F'(t) = f(t) then we can divide through to get a constant coefficient system which might be homogeneous or not. This is possibly preferable to using the type3/type4 solvers even though equation (2) will have a commuting antiderivative. If it is non-homogeneous then we will need to be able to replace t with tau which requires the inverse of F.

The 2nd order version of (2) is

X'' = f(t) A X + b(t)       (3)

Again we use X(t) = Y(tau) = Y(F(t)) to get

X'' = d/dt X' = d/dt (Y' F'(t)) = Y' F''(t) + Y'' F'(t)^2

Now the transformed equation is

F'(t)^2 Y'' + F''(t) Y' = f(t) A Y + b(t)

This might simplify for some cases of f and F. The old solvers have type6 and type10 for systems like (3). The type6 solver claims to be able to solve any system of this form but according to the tests it is broken and doesn't actually work.

The type10 solver solves a particular case where b(t)=0 and f(t) = 1 / (at^2 + bt + c)^2 with the substitution tau = F(t) = integral(1/(at^2 + bt + c), t). In that case F'(t) = 1 / (at^2 + bt + c) and F'(t)^2 = f(t). Also

F''(t) = d/dt (1 / (at^2 + bt + c)) = - (2at + b) / (at^2 + bt + c)^2 = -(2at + b) f(t)

So the transformed system would be

Y'' - (2at+b) Y' = A Y

This would be constant coefficient only if a is zero.

Instead the method suggested in the type10 solver is the substitution tau = F(t) and X(t) = Y(tau) f(t)^(1/4) but I'm not sure how that leads to a constant coefficient system as claimed. I can't see any tests for that solver.

Are there any working tests or examples for either type6 or type10?

[oscarbenjamin]

I think that the type6 solver can't possibly work in general. The class of system (3) as a non-homogeneous system would always need the substitution used for the type10 solver and I don't see how that can work except where f(t) is the inverse square of a quadratic i.e. f(t) = 1 / (a*t**2 + b*t + c).

In [258]: X = Function('X')                                                                                                                    

In [259]: Y, F, f = symbols('Y, F, f', cls=Function)                                                                                           

In [260]: g = Function('g')                                                                                                                    

In [261]: A = Symbol('A')                                                                                                                      

In [262]: eq = Eq(X(t).diff(t, 2), f(t) * A * X(t))                                                                                            

This is the system we want to solve:

In [263]: eq                                                                                                                                   
Out[263]: 
  2                    
 d                     
───(X(t)) = A⋅X(t)⋅f(t)
  2                    
dt                     

We'll take a general substitution of the form tau = F(t) and X(t) = Y(tau) * g(t):

In [264]: eq2 = eq.subs(X(t), Y(F(t)) * g(t)).doit()                                                                                           

In [265]: eq2                                                                                                                                  
Out[265]: 
⎛          2    2                2                     ⎞                  2                                                                 
⎜⎛d       ⎞    d                d          d           ⎟                 d            d          d            d                             
⎜⎜──(F(t))⎟ ⋅──────(Y(F(t))) + ───(F(t))⋅─────(Y(F(t)))⎟⋅g(t) + Y(F(t))⋅───(g(t)) + 2⋅──(F(t))⋅─────(Y(F(t)))⋅──(g(t)) = A⋅Y(F(t))⋅f(t)⋅g(t)
⎜⎝dt      ⎠       2              2       dF(t)         ⎟                  2           dt       dF(t)          dt                            
⎝            dF(t)             dt                      ⎠                dt                                                                  

In [266]: eq3 = Eq(*[s / (f(t)*g(t)) for s in eq2.args])                                                                                       

In [267]: eq3                                                                                                                                  
Out[267]: 
⎛          2    2                2                     ⎞                  2                                                       
⎜⎛d       ⎞    d                d          d           ⎟                 d            d          d            d                   
⎜⎜──(F(t))⎟ ⋅──────(Y(F(t))) + ───(F(t))⋅─────(Y(F(t)))⎟⋅g(t) + Y(F(t))⋅───(g(t)) + 2⋅──(F(t))⋅─────(Y(F(t)))⋅──(g(t))            
⎜⎝dt      ⎠       2              2       dF(t)         ⎟                  2           dt       dF(t)          dt                  
⎝            dF(t)             dt                      ⎠                dt                                                        
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── = A⋅Y(F(t))
                                                      f(t)⋅g(t)                                                                   

We now have a constant coeficient rhs so we need to make the coefficients on the lhs constant and we need the coefficient of Y' to be zero. First take the coefficient of Y'':

In [268]: lhs = eq3.lhs.expand().collect(Y(F(t)))                                                                                              

In [269]: lhs                                                                                                                                  
Out[269]: 
⎛  2                            ⎞                            2                   2    2           
⎜ d                             ⎟                           d          ⎛d       ⎞    d            
⎜───(F(t))     d        d       ⎟                  Y(F(t))⋅───(g(t))   ⎜──(F(t))⎟ ⋅──────(Y(F(t)))
⎜  2         2⋅──(F(t))⋅──(g(t))⎟                            2         ⎝dt      ⎠       2         
⎜dt            dt       dt      ⎟   d                      dt                      dF(t)          
⎜───────── + ───────────────────⎟⋅─────(Y(F(t))) + ───────────────── + ───────────────────────────
⎝   f(t)          f(t)⋅g(t)     ⎠ dF(t)                f(t)⋅g(t)                   f(t)           

In [270]: c2 = lhs.coeff(Y(F(t)).diff(F(t), 2))                                                                                                

In [271]: c2                                                                                                                                   
Out[271]: 
          2
⎛d       ⎞ 
⎜──(F(t))⎟ 
⎝dt      ⎠ 
───────────
    f(t)   

Given f(t) as specified by the input problem then we can choose F(t) for the transformation to make this coefficient 1 which gives

In [272]: dsolve(Eq(c2, 1), F(t))                                                                                                              
Out[272]: 
⎡            ⌠                          ⌠            ⎤
⎢            ⎮   ______                 ⎮   ______   ⎥
⎢F(t) = C₁ - ⎮ ╲╱ f(t)  dt, F(t) = C₁ + ⎮ ╲╱ f(t)  dt⎥
⎣            ⌡                          ⌡            ⎦

Using F(t) = Integral(sqrt(f(t)), t) we have F'(t) = sqrt(f(t)) the coefficient of Y' becomes:

In [275]: c1 = lhs.subs(F(t).diff(t), sqrt(f(t))).coeff(Y(F(t)).diff(F(t)))                                                                    

In [276]: c1                                                                                                                                   
Out[276]: 
d ⎛  ______⎞       d        
──⎝╲╱ f(t) ⎠     2⋅──(g(t)) 
dt                 dt       
──────────── + ─────────────
    f(t)         ______     
               ╲╱ f(t) ⋅g(t)

We can choose g(t) to make this zero:

In [277]: dsolve(c1, g(t))                                                                                                                     
Out[277]: 
          C₁   
g(t) = ────────
       4 ______
       ╲╱ f(t) 

We've now had to choose both F(t) and g(t) for the transformation and the coefficients of Y'' and Y' are 1 and 0. The coefficient of Y needs to be constant:

In [281]: rep = {F(t).diff(t):sqrt(f(t)), g(t):1/f(t)**Rational(1, 4)}                                                                         

In [282]: lhs.subs(rep).doit()                                                                                                                 
Out[282]: 
⎛                            2⎞                          
⎜                  ⎛d       ⎞ ⎟                          
⎜      2         5⋅⎜──(f(t))⎟ ⎟                          
⎜     d            ⎝dt      ⎠ ⎟                          
⎜- 4⋅───(f(t)) + ─────────────⎟⋅Y(F(t))                  
⎜      2              f(t)    ⎟              2           
⎝    dt                       ⎠             d            
─────────────────────────────────────── + ──────(Y(F(t)))
                    2                          2         
                16⋅f (t)                  dF(t)          

In [283]: c0 = lhs.subs(rep).doit().coeff(Y(F(t)))                                                                                             

In [284]: c0                                                                                                                                   
Out[284]: 
                            2
                  ⎛d       ⎞ 
      2         5⋅⎜──(f(t))⎟ 
     d            ⎝dt      ⎠ 
- 4⋅───(f(t)) + ─────────────
      2              f(t)    
    dt                       
─────────────────────────────
               2             
           16⋅f (t)          

This coefficient will be constant as given by the particular case for f(t) specified for the type10 solver:

In [285]: a, b, c = symbols('a, b, c')                                                                                                         

In [286]: c0.subs(f(t), 1/(a*t**2 + b*t + c)**2).doit().simplify()                                                                             
Out[286]: 
       2
      b 
a⋅c - ──
      4 

Other choices of simple functions do not lead to this being constant:

In [294]: c0.diff(t).subs(f(t), t).doit().simplify()                                                                                           
Out[294]: 
 -15 
─────
    4
16⋅t 

In [295]: c0.diff(t).subs(f(t), exp(t)).doit().simplify()                                                                                      
Out[295]: 
  -t 
-ℯ   
─────
  16 

We can view c0.diff(t) as an ODE whose solutions are the possible choices for f(t) such that we can use this kind of substitution to get to a constant coefficient 2nd order ODE. It's a 3rd order nonlinear ODE:

In [301]: c0.diff(t).as_numer_denom()[0].expand()                                                                                              
Out[301]: 
            3                             2                            3
     7     d              6    d         d              5    ⎛d       ⎞ 
- 4⋅f (t)⋅───(f(t)) + 18⋅f (t)⋅──(f(t))⋅───(f(t)) - 15⋅f (t)⋅⎜──(f(t))⎟ 
            3                  dt         2                  ⎝dt      ⎠ 
          dt                            dt                              

I don't know how we could go about getting an analytic solution for this ODE (SymPy can't and Wolfram Alpha can't). We can see that f(t) = 1 / (a*t**2 + b*t + c) is a solution of this by substituting it in though. In fact being a 3-parameter family of solutions to a 3rd order ODE it possibly represents the general solution.

If indeed 1 / (a*t**2 + b*t + c) is the general solution of Out[301] then I think that means this class of substitution can not work for a system that has any other kind of function f(t).

In that case I think it is reasonable to implement the type10 solver is a particular substitution to be done with systems of 2nd order ODEs. It can however be generalised to arbitrary numbers of equations rather than just a pair of 2nd order ODEs.

[mijo2]

In that case I think it is reasonable to implement the type10 solver is a particular substitution to be done with systems of 2nd order ODEs. It can however be generalised to arbitrary numbers of equations rather than just a pair of 2nd order ODEs.

Well it looks as if the function f(t) = 1 / (a*t**2 + b*t + c) maybe is the solution for this ode. So, for now, I would implement the method which will simply extend the work done in type 10 and add a note referring to maybe this comment or issue so that if someone wants to extend this method, they will have a good starting point.

Thanks for all your analysis. I would implement the following for now when it comes to higher-order systems:

1. Complete the implementation of the Cauchy-Euler technique
2. Add the method for type 10 substitution
3. Add appropriate test cases and find issues solved by our methods
4. Add a specific method wherein the system will take the information of the minimum and maximum orders of the dependent variables and if, one of the dependent variables has an order like: max(x(t)) > min(x(t)) > 1, then we can safely reduce the system beforehand, so that we don't have unnecessarily big matrices.

[oscarbenjamin]

Well it looks as if the function f(t) = 1 / (a*t**2 + b*t + c) maybe is the solution for this ode.

Yes actually we can show it by substituting f(t) = 1/g(t)**2:

In [106]: eq = -4 * f(t)**7 * f(t).diff(t, 3) + 18 * f(t)**6 * f(t).diff(t) * f(t).diff(t, 2) - 15 * f(t)**5 * (f(t).diff(t))**3               

In [107]: eq                                                                                                                                   
Out[107]: 
            3                             2                            3
     7     d              6    d         d              5    ⎛d       ⎞ 
- 4⋅f (t)⋅───(f(t)) + 18⋅f (t)⋅──(f(t))⋅───(f(t)) - 15⋅f (t)⋅⎜──(f(t))⎟ 
            3                  dt         2                  ⎝dt      ⎠ 
          dt                            dt                              

In [108]: eq.subs(f(t), 1/g(t)**2)                                                                                                             
Out[108]: 
                3                  2              3       
     ⎛d ⎛  1  ⎞⎞       d ⎛  1  ⎞  d ⎛  1  ⎞      d ⎛  1  ⎞
  15⋅⎜──⎜─────⎟⎟    18⋅──⎜─────⎟⋅───⎜─────⎟   4⋅───⎜─────⎟
     ⎜dt⎜ 2   ⎟⎟       dt⎜ 2   ⎟   2⎜ 2   ⎟       3⎜ 2   ⎟
     ⎝  ⎝g (t)⎠⎠         ⎝g (t)⎠ dt ⎝g (t)⎠     dt ⎝g (t)⎠
- ─────────────── + ─────────────────────── - ────────────
        10                    12                  14      
       g  (t)                g  (t)              g  (t)   

In [109]: eq.subs(f(t), 1/g(t)**2).doit().expand()                                                                                             
Out[109]: 
    3      
   d       
8⋅───(g(t))
    3      
  dt       
───────────
    17     
   g  (t)  

In [110]: dsolve(_)                                                                                                                            
Out[110]: 
                       2
g(t) = C₁ + C₂⋅t + C₃⋅t 

In [111]: 1/_.rhs**2                                                                                                                           
Out[111]: 
         1          
────────────────────
                   2
⎛                2⎞ 
⎝C₁ + C₂⋅t + C₃⋅t ⎠ 

I would implement the following

Yes, that sounds good.

[mijo2]

If it is non-homogeneous then we will need to be able to replace t with tau which requires the inverse of F.

If tau doesn't have an inverse, then probably we can revert back to the old method of solving.

[mijo2]

I think that the type6 solver can't possibly work in general.

This is true, but it looks like its probably able to solve a system of the form:

X'' = A*X

For example(tested on master):

In [12]: eq8                                                                                                                                                                                                
Out[12]: 
(Eq(Derivative(x(t), (t, 2)), t*(4*x(t) + 9*y(t))),
 Eq(Derivative(y(t), (t, 2)), t*(12*x(t) - 6*y(t))))

In [13]: sol = dsolve(eq8)                                                                                                                                                                                  

In [14]: checksysodesol(eq8, sol)                                                                                                                                                                           
Out[14]: (True, [0, 0])

In [15]: classify_sysode(eq8)                                                                                                                                                                               
Out[15]: 
{'no_of_equation': 2,
 'eq': [-t*(4*x(t) + 9*y(t)) + Derivative(x(t), (t, 2)),
  -t*(12*x(t) - 6*y(t)) + Derivative(y(t), (t, 2))],
 'func': [x(t), y(t)],
 'order': {x(t): 2, y(t): 2},
 'func_coeff': {(0, x(t), 0): -4*t,
  (0, x(t), 1): 0,
  (0, x(t), 2): 1,
  (0, y(t), 0): -9*t,
  (0, y(t), 1): 0,
  (0, y(t), 2): 0,
  (1, x(t), 0): -12*t,
  (1, x(t), 1): 0,
  (1, x(t), 2): 0,
  (1, y(t), 0): 6*t,
  (1, y(t), 1): 0,
  (1, y(t), 2): 1},
 'is_linear': True,
 'type_of_equation': 'type6'}

Now for this particular form, you have already mentioned a technique to be implemented, but what if there are other types of systems that type6 solver can solve, not all but what if by removing this solver, we may not have any other solver to handle those cases for 2 equations, Order 2. How can we identify for which f(t), type 6 solver actually works? I would try to find a method or maybe manually try some trivial systems.

[oscarbenjamin]

Okay, I guess I didn't look closely at the type6 solver. It's using another technique based on diagonalisation I think. We don't know that the matrix is always diagonalisable so let's use the Jordan normal form. We have a system like

X'' = f(t) M X

Suppose that the JNF of M is P J P^-1 then

X'' = f(t) P J P^-1 X
P^-1 X'' = f(t) J P^-1 X

Now since M is constant P and J are as well. We can use the substitution Z = P^-1 X with Z'' = P^-1 X'' which gives

Z'' = f(t) J Z

Now if M was diagonalisable then J is diagonal and the system for Z is completely decoupled: the weakly connected components are single ODEs. If M was not diagonalisable then the weakly connected components are the Jordan blocks and the strongly connected components of that are the individual equations for each z_i. So we can solve this as a bunch of individual ODEs for Z and then we find X as X = P Z.

Here's how this works for the example you gave:

In [16]: eq = (Eq(Derivative(x(t), (t, 2)), t*(4*x(t) + 9*y(t))), 
    ...:  Eq(Derivative(y(t), (t, 2)), t*(12*x(t) - 6*y(t))))                                                                                  

In [17]: eq                                                                                                                                    
Out[17]: 
⎛  2                                2                             ⎞
⎜ d                                d                              ⎟
⎜───(x(t)) = t⋅(4⋅x(t) + 9⋅y(t)), ───(y(t)) = t⋅(12⋅x(t) - 6⋅y(t))⎟
⎜  2                                2                             ⎟
⎝dt                               dt                              ⎠

In [18]: from sympy.solvers.ode import linear_ode_to_matrix                                                                                    

In [19]: linear_ode_to_matrix(eq, [x(t), y(t)], t, 2)                                                                                          
Out[19]: 
⎛⎡⎡1  0⎤  ⎡0  0⎤  ⎡4⋅t   9⋅t ⎤⎤  ⎡0⎤⎞
⎜⎢⎢    ⎥, ⎢    ⎥, ⎢          ⎥⎥, ⎢ ⎥⎟
⎝⎣⎣0  1⎦  ⎣0  0⎦  ⎣12⋅t  -6⋅t⎦⎦  ⎣0⎦⎠

In [20]: M = linear_ode_to_matrix(eq, [x(t), y(t)], t, 2)[0][2]                                                                                

In [21]: Mp = M / t                                                                                                                            

In [22]: Mp                                                                                                                                    
Out[22]: 
Matrix([
[ 4,  9],
[12, -6]])

In [23]: P, J = Mp.jordan_form()                                                                                                               

In [24]: P                                                                                                                                     
Out[24]: 
Matrix([
[-9/(5 - sqrt(133)), -9/(5 + sqrt(133))],
[                 1,                  1]])

In [25]: J                                                                                                                                     
Out[25]: 
Matrix([
[-1 + sqrt(133),              0],
[             0, -sqrt(133) - 1]])

In [26]: z1, z2 = symbols('z1, z2', cls=Function)                                                                                              

In [27]: Z = Matrix([z1(t), z2(t)])                                                                                                            

In [28]: t*J*Z - Z.diff(t, 2)                                                                                                                  
Out[28]: 
Matrix([
[t*(-1 + sqrt(133))*z1(t) - Derivative(z1(t), (t, 2))],
[t*(-sqrt(133) - 1)*z2(t) - Derivative(z2(t), (t, 2))]])

In [29]: dsolve(t*J*Z - Z.diff(t, 2))                                                                                                          
---------------------------------------------------------------------------
NotImplementedError

I tested that with this PR. I think that there is maybe a problem with the component solvers. Do they pass single ODEs back to dsolve? That example should be solvable by splitting into weakly connected components.

We can solve it manually:

In [33]: eqz1, eqz2 = t*J*Z - Z.diff(t, 2)                                                                                                     

In [34]: sol1 = dsolve(eqz1)                                                                                                                   

In [35]: sol1                                                                                                                                  
Out[35]: 
             ⎛  3 ___________⎞        ⎛  3 ___________⎞
z₁(t) = C₁⋅Ai⎝t⋅╲╱ -1 + √133 ⎠ + C₂⋅Bi⎝t⋅╲╱ -1 + √133 ⎠

In [36]: sol2 = dsolve(eqz2)                                                                                                                   

In [37]: sol2                                                                                                                                  
Out[37]: 
             ⎛   3 __________⎞        ⎛   3 __________⎞
z₂(t) = C₁⋅Ai⎝-t⋅╲╱ 1 + √133 ⎠ + C₂⋅Bi⎝-t⋅╲╱ 1 + √133 ⎠

In [39]: sol2 = sol2.subs(zip(symbols('C1:3'), symbols('C3:5')))                                                                               

In [40]: sol2                                                                                                                                  
Out[40]: 
             ⎛   3 __________⎞        ⎛   3 __________⎞
z₂(t) = C₃⋅Ai⎝-t⋅╲╱ 1 + √133 ⎠ + C₄⋅Bi⎝-t⋅╲╱ 1 + √133 ⎠

In [41]: solZ = Matrix([sol1.rhs, sol2.rhs])                                                                                                   

In [42]: solZ                                                                                                                                  
Out[42]: 
Matrix([
[C1*airyai(t*(-1 + sqrt(133))**(1/3)) + C2*airybi(t*(-1 + sqrt(133))**(1/3))],
[C3*airyai(-t*(1 + sqrt(133))**(1/3)) + C4*airybi(-t*(1 + sqrt(133))**(1/3))]])

In [43]: P * solZ                                                                                                                              
Out[43]: 
Matrix([
[-9*(C1*airyai(t*(-1 + sqrt(133))**(1/3)) + C2*airybi(t*(-1 + sqrt(133))**(1/3)))/(5 - sqrt(133)) - 9*(C3*airyai(-t*(1 + sqrt(133))**(1/3)) + C4*airybi(-t*(1 + sqrt(133))**(1/3)))/(5 + sqrt(133))],
[                                         C1*airyai(t*(-1 + sqrt(133))**(1/3)) + C2*airybi(t*(-1 + sqrt(133))**(1/3)) + C3*airyai(-t*(1 + sqrt(133))**(1/3)) + C4*airybi(-t*(1 + sqrt(133))**(1/3))]])

In [44]: solX = list(P * solZ)                                                                                                                 

In [45]: solX                                                                                                                                  
Out[45]: 
⎡    ⎛     ⎛  3 ___________⎞        ⎛  3 ___________⎞⎞     ⎛     ⎛   3 __________⎞        ⎛   3 __________⎞⎞                                 
⎢  9⋅⎝C₁⋅Ai⎝t⋅╲╱ -1 + √133 ⎠ + C₂⋅Bi⎝t⋅╲╱ -1 + √133 ⎠⎠   9⋅⎝C₃⋅Ai⎝-t⋅╲╱ 1 + √133 ⎠ + C₄⋅Bi⎝-t⋅╲╱ 1 + √133 ⎠⎠       ⎛  3 ___________⎞        ⎛
⎢- ─────────────────────────────────────────────────── - ───────────────────────────────────────────────────, C₁⋅Ai⎝t⋅╲╱ -1 + √133 ⎠ + C₂⋅Bi⎝
⎣                        5 - √133                                              5 + √133                                                      

                                                                  ⎤
  3 ___________⎞        ⎛   3 __________⎞        ⎛   3 __________⎞⎥
t⋅╲╱ -1 + √133 ⎠ + C₃⋅Ai⎝-t⋅╲╱ 1 + √133 ⎠ + C₄⋅Bi⎝-t⋅╲╱ 1 + √133 ⎠⎥
                                                                  ⎦

In [46]: solX = [Eq(x(t), solX[0]), Eq(y(t), solX[1])]                                                                                         

In [47]: solX                                                                                                                                  
Out[47]: 
⎡           ⎛     ⎛  3 ___________⎞        ⎛  3 ___________⎞⎞     ⎛     ⎛   3 __________⎞        ⎛   3 __________⎞⎞                          
⎢         9⋅⎝C₁⋅Ai⎝t⋅╲╱ -1 + √133 ⎠ + C₂⋅Bi⎝t⋅╲╱ -1 + √133 ⎠⎠   9⋅⎝C₃⋅Ai⎝-t⋅╲╱ 1 + √133 ⎠ + C₄⋅Bi⎝-t⋅╲╱ 1 + √133 ⎠⎠              ⎛  3 _______
⎢x(t) = - ─────────────────────────────────────────────────── - ───────────────────────────────────────────────────, y(t) = C₁⋅Ai⎝t⋅╲╱ -1 + √
⎣                               5 - √133                                              5 + √133                                               

                                                                                ⎤
____⎞        ⎛  3 ___________⎞        ⎛   3 __________⎞        ⎛   3 __________⎞⎥
133 ⎠ + C₂⋅Bi⎝t⋅╲╱ -1 + √133 ⎠ + C₃⋅Ai⎝-t⋅╲╱ 1 + √133 ⎠ + C₄⋅Bi⎝-t⋅╲╱ 1 + √133 ⎠⎥
                                                                                ⎦

In [48]: checkodesol(eq, solX)                                                                                                                 
Out[48]: (True, [0, 0])

[mijo2]

Do they pass single ODEs back to dsolve? That example should be solvable by splitting into weakly connected components.

No, they don't pass the single ODEs back to dsolve. The weakly connected components part work fine, but the matrix of a second order system with only one dependent variable and which is of the form: X'' = AX is not commutative with its antiderivative if reduced first order by introducing dummy functions. If we add a dsolve call for single equation systems, then that may help extend dsolve_system's capabilities.

[mijo2]

Also 1st order systems of the form

X' = f(t) A X + b(t) (2)

Should I introduce new types, like type5 and type6 for this particular case, since we are using a different solution?

Also, it would be for the best if we use this solution only when F(t) = integrate(f(t), t) is invertible. I tried a solution when F(t) = t**2/2, but it gave incorrect results for both of the substitutions of its inverses.

[oscarbenjamin]

Should I introduce new types, like type5 and type6 for this particular case

You could do. We need to decide whether it's worth having them as special cases though.

Any system of the form (2) will have a commuting antiderivative so it can already be solved. The method described above using JNF for 2nd order systems also works for 1st (or any) order so again the question is really which gives the best (simplest, fastest) result.

[mijo2]

@oscarbenjamin I tested with some invertible functions(using solveset) and found that most of the times, the new solution is approximately 2-3 times faster than the old "type3" and "type4" solvers

[oscarbenjamin]

2-3 times faster

Does it give the solution in the same form or is it better in some way?

I think that the most important thing to finish before the end of GSOC is simplification of the solutions from the new solvers. Have you made any progress on that?

Adding solvers for more types of system can be done later as it is something that can be done incrementally. The simplification of the solutions is something that needs to be resolved before the next release of sympy. I think there are some cases that used to be handled by dsolve and still are but now have more complicated solutions.

Integrals should be evaluated by default although there should be an option to disable that. Exponential factors should be combined like exp(t*(1 + sqrt(2))) rather than exp(t)*exp(t*sqrt(2)) or exp(t + t*sqrt(2)). Also we need to collect on factors involving the dependent variable and on integration constants.

For example sympy 1.6 gives:

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

In [2]: eq = [Eq(x(t).diff(t), x(t) + y(t) + 1), Eq(y(t).diff(t), x(t) + y(t))]                                                                               

In [3]: dsolve(eq)                                                                                                                                            
Out[3]: 
⎡           2⋅t        t   1             2⋅t        t   1⎤
⎢x(t) = C₁⋅ℯ    + C₂ + ─ - ─, y(t) = C₁⋅ℯ    - C₂ - ─ - ─⎥
⎣                      2   4                        2   4⎦

However on master we have:

In [4]: dsolve(eq)                                                                                                                                            
Out[4]: 
⎡                            ⌠                                                  ⌠                     ⎤
⎢                            ⎮  -2⋅t                                            ⎮  -2⋅t               ⎥
⎢                 2⋅t    2⋅t ⎮ ℯ          ⌠                          2⋅t    2⋅t ⎮ ℯ          ⌠        ⎥
⎢x(t) = -C₁ + C₂⋅ℯ    + ℯ   ⋅⎮ ───── dt - ⎮ -1/2 dt, y(t) = C₁ + C₂⋅ℯ    + ℯ   ⋅⎮ ───── dt + ⎮ -1/2 dt⎥
⎢                            ⎮   2        ⌡                                     ⎮   2        ⌡        ⎥
⎣                            ⌡                                                  ⌡                     ⎦

That's just a case of calling doit. Some of the more complicated examples come from the type3/4 solvers:

In [9]: eq = [Eq(x(t).diff(t), t*x(t) + t*y(t) + 1), Eq(y(t).diff(t), t*x(t) + t*y(t))]                                                                       

In [10]: dsolve(eq)                                                                                                                                           
Out[10]: 
⎡                                             ⌠                         ⌠                 ⌠                   ⌠                                             
⎢                                             ⎮ ⎛         2⎞            ⎮ ⎛       2⎞      ⎮ ⎛         2⎞      ⎮ ⎛       2⎞                                  
⎢                                        ⎛ 2⎞ ⎮ ⎜       -t ⎟       ⎛ 2⎞ ⎮ ⎜     -t ⎟      ⎮ ⎜       -t ⎟      ⎮ ⎜     -t ⎟                                  
⎢                                        ⎝t ⎠ ⎮ ⎜  1   ℯ   ⎟       ⎝t ⎠ ⎮ ⎜1   ℯ   ⎟      ⎮ ⎜  1   ℯ   ⎟      ⎮ ⎜1   ℯ   ⎟                                  
⎢           ⎛ 2⎞            ⎛ 2⎞        ℯ    ⋅⎮ ⎜- ─ + ────⎟ dt   ℯ    ⋅⎮ ⎜─ + ────⎟ dt   ⎮ ⎜- ─ + ────⎟ dt   ⎮ ⎜─ + ────⎟ dt             ⎛ 2⎞            ⎛ 
⎢           ⎝t ⎠            ⎝t ⎠              ⎮ ⎝  2    2  ⎠            ⎮ ⎝2    2  ⎠      ⎮ ⎝  2    2  ⎠      ⎮ ⎝2    2  ⎠                ⎝t ⎠            ⎝t
⎢       C₁⋅ℯ       C₁   C₂⋅ℯ       C₂         ⌡                         ⌡                 ⌡                   ⌡                       C₁⋅ℯ       C₁   C₂⋅ℯ  
⎢x(t) = ──────── - ── + ──────── + ── + ─────────────────────── + ───────────────────── - ───────────────── + ───────────────, y(t) = ──────── + ── + ──────
⎣          2       2       2       2               2                        2                     2                  2                   2       2       2  

                ⌠                         ⌠                 ⌠                   ⌠              ⎤
                ⎮ ⎛         2⎞            ⎮ ⎛       2⎞      ⎮ ⎛         2⎞      ⎮ ⎛       2⎞   ⎥
           ⎛ 2⎞ ⎮ ⎜       -t ⎟       ⎛ 2⎞ ⎮ ⎜     -t ⎟      ⎮ ⎜       -t ⎟      ⎮ ⎜     -t ⎟   ⎥
           ⎝t ⎠ ⎮ ⎜  1   ℯ   ⎟       ⎝t ⎠ ⎮ ⎜1   ℯ   ⎟      ⎮ ⎜  1   ℯ   ⎟      ⎮ ⎜1   ℯ   ⎟   ⎥
2⎞        ℯ    ⋅⎮ ⎜- ─ + ────⎟ dt   ℯ    ⋅⎮ ⎜─ + ────⎟ dt   ⎮ ⎜- ─ + ────⎟ dt   ⎮ ⎜─ + ────⎟ dt⎥
 ⎠              ⎮ ⎝  2    2  ⎠            ⎮ ⎝2    2  ⎠      ⎮ ⎝  2    2  ⎠      ⎮ ⎝2    2  ⎠   ⎥
     C₂         ⌡                         ⌡                 ⌡                   ⌡              ⎥
── - ── + ─────────────────────── + ───────────────────── + ───────────────── - ───────────────⎥
     2               2                        2                     2                  2       ⎦

There are ways to simplify that after evaluating the integrals by using collect:

In [25]: solx, soly = dsolve(eq)                                                                                                                              

In [26]: solx = solx.doit()                                                                                                                                   

In [27]: solx                                                                                                                                                 
Out[27]: 
                                                              ⎛ 2⎞                    ⎛ 2⎞
           ⎛ 2⎞            ⎛ 2⎞            ⎛  t   √π⋅erf(t)⎞  ⎝t ⎠   ⎛t   √π⋅erf(t)⎞  ⎝t ⎠
           ⎝t ⎠            ⎝t ⎠            ⎜- ─ + ─────────⎟⋅ℯ       ⎜─ + ─────────⎟⋅ℯ    
       C₁⋅ℯ       C₁   C₂⋅ℯ       C₂   t   ⎝  2       4    ⎠         ⎝2       4    ⎠      
x(t) = ──────── - ── + ──────── + ── + ─ + ─────────────────────── + ─────────────────────
          2       2       2       2    2              2                        2          

In [28]: solx.rhs.collect(symbols('C1:3')).collect(solx.atoms(exp))                                                                                           
Out[28]: 
   ⎛ ⎛ 2⎞    ⎞      ⎛ ⎛ 2⎞    ⎞           ⎛ 2⎞       
   ⎜ ⎝t ⎠    ⎟      ⎜ ⎝t ⎠    ⎟           ⎝t ⎠       
   ⎜ℯ       1⎟      ⎜ℯ       1⎟   t   √π⋅ℯ    ⋅erf(t)
C₁⋅⎜───── - ─⎟ + C₂⋅⎜───── + ─⎟ + ─ + ───────────────
   ⎝  2     2⎠      ⎝  2     2⎠   2          4       

In [29]: solx.rhs.collect(solx.atoms(exp))                                                                                                                    
Out[29]: 
                                       ⎛ 2⎞
  C₁   C₂   t   ⎛C₁   C₂   √π⋅erf(t)⎞  ⎝t ⎠
- ── + ── + ─ + ⎜── + ── + ─────────⎟⋅ℯ    
  2    2    2   ⎝2    2        4    ⎠ 

I'm not sure which of those forms is best. I quite like seeing the constants brought out to the top in Out[28] as it shows the fundamental linearly independent solutions and the non-homogeneous term. On the other hand Out[29] makes the overall dependence on t a bit clearer. An important consideration is which form leads to the smallest expressions for the solution in the cases where the solution is big and complicated.

It's possible that powsimp, factor_terms and/or expand_mul can be useful in this as well.

We also need to do something about the large fragile tests that were created using solsimp/simpsol.

[mijo2]

Does it give the solution in the same form or is it better in some way?

Almost the similar form.

[oscarbenjamin]

I have made a function for simplifying the solutions of dsolve_system. It helps a lot in the simple cases but is variable for complicated cases:

def simpsol(sol, wrt1, wrt2, doit=True):
    """Simplify solutions from dsolve_system."""

    # The parameter sol is the solution as returned by dsolve (list of Eq).
    #
    # The parameters wrt1 and wrt2 are lists of symbols to be collected for
    # with those in wrt1 being collected for first. This allows for collecting
    # on any factors involving the independent variable before collecting on
    # the integration constants or vice versa using e.g.:
    #
    #     sol = simpsol(sol, [t], [C1, C2])  # t first, constants after
    #     sol = simpsol(sol, [C1, C2], [t])  # constants first, t after
    #
    # If doit=True (default) then simpsol will begin by evaluating any
    # unevaluated integrals. Since many integrals will appear multiple times
    # in the solutions this is done intelligently by computing each integral
    # only once.
    #
    # The strategy is to first perform simple cancellation with factor_terms
    # and then multiply out all brackets with expand_mul. This gives an Add
    # with many terms.
    #
    # We split each term into two multiplicative factors dep and coeff where
    # all factors that involve wrt1 are in dep and any constant factors are in
    # coeff e.g.
    #         sqrt(2)*C1*exp(t) -> ( exp(t) , sqrt(2)*C1 )
    #
    # The dep factors are simplified using powsimp to combine expanded
    # exponential factors e.g.
    #              exp(a*t)*exp(b*t) -> exp(t*(a+b))
    #
    # We then collect coefficients for all terms having the same (simplified)
    # dep. The coefficients are then simplified using together and ratsimp and
    # lastly by recursively applying the same transformation to the
    # coefficients to collect on wrt2.
    #
    # Finally the result is recombined into an Add and signsimp is used to
    # normalise any minus signs.

    def simprhs(rhs, rep, wrt1, wrt2):
        """Simplify the rhs of an ODE solution"""
        if rep:
            rhs = rhs.subs(rep)
        rhs = factor_terms(rhs)
        rhs = simp_coeff_dep(rhs, wrt1, wrt2)
        rhs = signsimp(rhs)
        return rhs

    def simp_coeff_dep(expr, wrt1, wrt2=None):
        """Split rhs into terms, split terms into dep and coeff and collect on dep"""
        terms = Add.make_args(expand_mul(expr))
        dc = {}
        for term in terms:
            coeff, dep = term.as_independent(*wrt1, as_Add=False)
            # Collect together the coefficients for terms that have the same
            # dependence on wrt1 (after dep is normalised using simpdep).
            dep = simpdep(dep, wrt1)
            if dep not in dc:
                dc[dep] = coeff
            else:
                dc[dep] += coeff
        # Apply the method recursively to the coefficients but this time
        # collecting on wrt2 rather than wrt2.
        termpairs = ((simpcoeff(c, wrt2), d) for d, c in dc.items())
        if wrt2 is not None:
            termpairs = ((simp_coeff_dep(c, wrt2), d) for c, d in termpairs)
        return Add(*(c * d for c, d in termpairs))

    def simpdep(term, wrt1):
        """Normalise factors involving t with powsimp and recombine exp"""
        def canonicalise(a):
            # Using factor_terms here isn't quite right because it leads to things
            # like exp(t*(1+t)) that we don't want. We do want to cancel factors
            # and pull out a common denominator but ideally the numerator would be
            # expressed as a standard form polynomial in t so we expand_mul
            # and collect afterwards.
            a = factor_terms(a)
            num, den = a.as_numer_denom()
            num = expand_mul(num)
            num = collect(num, wrt1)
            if den != 1:
                return Mul(num, Pow(den, -1, evaluate=False), evaluate=False)
            else:
                return num

        term = powsimp(term)
        rep = {e: exp(canonicalise(e.args[0])) for e in term.atoms(exp)}
        term = term.subs(rep)
        return term

    def simpcoeff(coeff, wrt2):
        """Bring to a common fraction and cancel with ratsimp"""
        coeff = together(coeff)
        if coeff.is_polynomial():
            # Calling ratsimp can be expensive. The main reason is to simplify
            # sums of terms with irrational denominators so we limit ourselves
            # to the case where the expression is polynomial in any symbols.
            # Maybe there's a better approach...
            coeff = ratsimp(radsimp(coeff))
        # collect on secondary variables first and any remaining symbols after
        if wrt2 is not None:
            syms = list(wrt2) + list(ordered(coeff.free_symbols - set(wrt2)))
        else:
            syms = list(ordered(coeff.free_symbols))
        coeff = collect(coeff, syms)
        coeff = together(coeff)
        return coeff

    # There are often repeated integrals. Collect unique integrals and
    # evaluate each once and then substitute into the final result to replace
    # all occurrences in each of the solution equations.
    if doit:
        integrals = set().union(*(s.atoms(Integral) for s in sol))
        rep = {i: factor_terms(i).doit() for i in integrals}
    else:
        rep = {}

    sol = [Eq(s.lhs, simprhs(s.rhs, rep, wrt1, wrt2)) for s in sol]

    return sol

That gives e.g.:

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

In [3]: eq = [Eq(x(t).diff(t), x(t) + y(t) + 1), Eq(y(t).diff(t), x(t) + y(t))]                                                                               

In [4]: sol = dsolve(eq)                                                                                                                                      

In [5]: sol                                                                                                                                                   
Out[5]: 
⎡                            ⌠                                                  ⌠                     ⎤
⎢                            ⎮  -2⋅t                                            ⎮  -2⋅t               ⎥
⎢                 2⋅t    2⋅t ⎮ ℯ          ⌠                          2⋅t    2⋅t ⎮ ℯ          ⌠        ⎥
⎢x(t) = -C₁ + C₂⋅ℯ    + ℯ   ⋅⎮ ───── dt - ⎮ -1/2 dt, y(t) = C₁ + C₂⋅ℯ    + ℯ   ⋅⎮ ───── dt + ⎮ -1/2 dt⎥
⎢                            ⎮   2        ⌡                                     ⎮   2        ⌡        ⎥
⎣                            ⌡                                                  ⌡                     ⎦

In [7]: simpsol(sol, [t], symbols('C1:3'))                                                                                                                    
Out[7]: 
⎡                 2⋅t   t   1                  2⋅t   t   1⎤
⎢x(t) = -C₁ + C₂⋅ℯ    + ─ - ─, y(t) = C₁ + C₂⋅ℯ    - ─ - ─⎥
⎣                       2   4                        2   4⎦

In [8]: eq = [Eq(x(t).diff(t), t*x(t) + t*y(t) + 1), Eq(y(t).diff(t), t*x(t) + t*y(t))]                                                                       

In [9]: sol = dsolve(eq)                                                                                                                                      

In [10]: sol                                                                                                                                                  
Out[10]: 
⎡                                             ⌠                         ⌠                 ⌠                   ⌠                                             
⎢                                             ⎮ ⎛         2⎞            ⎮ ⎛       2⎞      ⎮ ⎛         2⎞      ⎮ ⎛       2⎞                                  
⎢                                        ⎛ 2⎞ ⎮ ⎜       -t ⎟       ⎛ 2⎞ ⎮ ⎜     -t ⎟      ⎮ ⎜       -t ⎟      ⎮ ⎜     -t ⎟                                  
⎢                                        ⎝t ⎠ ⎮ ⎜  1   ℯ   ⎟       ⎝t ⎠ ⎮ ⎜1   ℯ   ⎟      ⎮ ⎜  1   ℯ   ⎟      ⎮ ⎜1   ℯ   ⎟                                  
⎢           ⎛ 2⎞            ⎛ 2⎞        ℯ    ⋅⎮ ⎜- ─ + ────⎟ dt   ℯ    ⋅⎮ ⎜─ + ────⎟ dt   ⎮ ⎜- ─ + ────⎟ dt   ⎮ ⎜─ + ────⎟ dt             ⎛ 2⎞            ⎛ 
⎢           ⎝t ⎠            ⎝t ⎠              ⎮ ⎝  2    2  ⎠            ⎮ ⎝2    2  ⎠      ⎮ ⎝  2    2  ⎠      ⎮ ⎝2    2  ⎠                ⎝t ⎠            ⎝t
⎢       C₁⋅ℯ       C₁   C₂⋅ℯ       C₂         ⌡                         ⌡                 ⌡                   ⌡                       C₁⋅ℯ       C₁   C₂⋅ℯ  
⎢x(t) = ──────── - ── + ──────── + ── + ─────────────────────── + ───────────────────── - ───────────────── + ───────────────, y(t) = ──────── + ── + ──────
⎣          2       2       2       2               2                        2                     2                  2                   2       2       2  

                ⌠                         ⌠                 ⌠                   ⌠              ⎤
                ⎮ ⎛         2⎞            ⎮ ⎛       2⎞      ⎮ ⎛         2⎞      ⎮ ⎛       2⎞   ⎥
           ⎛ 2⎞ ⎮ ⎜       -t ⎟       ⎛ 2⎞ ⎮ ⎜     -t ⎟      ⎮ ⎜       -t ⎟      ⎮ ⎜     -t ⎟   ⎥
           ⎝t ⎠ ⎮ ⎜  1   ℯ   ⎟       ⎝t ⎠ ⎮ ⎜1   ℯ   ⎟      ⎮ ⎜  1   ℯ   ⎟      ⎮ ⎜1   ℯ   ⎟   ⎥
2⎞        ℯ    ⋅⎮ ⎜- ─ + ────⎟ dt   ℯ    ⋅⎮ ⎜─ + ────⎟ dt   ⎮ ⎜- ─ + ────⎟ dt   ⎮ ⎜─ + ────⎟ dt⎥
 ⎠              ⎮ ⎝  2    2  ⎠            ⎮ ⎝2    2  ⎠      ⎮ ⎝  2    2  ⎠      ⎮ ⎝2    2  ⎠   ⎥
     C₂         ⌡                         ⌡                 ⌡                   ⌡              ⎥
── - ── + ─────────────────────── + ───────────────────── + ───────────────── - ───────────────⎥
     2               2                        2                     2                  2       ⎦

In [11]: simpsol(sol, [t], symbols('C1:3'))                                                                                                                   
Out[11]: 
⎡                                             ⎛ 2⎞                                                    ⎛ 2⎞       ⎤
⎢                                  ⎛ 2⎞       ⎝t ⎠                                         ⎛ 2⎞       ⎝t ⎠       ⎥
⎢         C₁   C₂   t   ⎛C₁   C₂⎞  ⎝t ⎠   √π⋅ℯ    ⋅erf(t)         C₁   C₂   t   ⎛C₁   C₂⎞  ⎝t ⎠   √π⋅ℯ    ⋅erf(t)⎥
⎢x(t) = - ── + ── + ─ + ⎜── + ──⎟⋅ℯ     + ───────────────, y(t) = ── - ── - ─ + ⎜── + ──⎟⋅ℯ     + ───────────────⎥
⎣         2    2    2   ⎝2    2 ⎠                4                2    2    2   ⎝2    2 ⎠                4       ⎦

In [12]: simpsol(sol, symbols('C1:3'), [t])                                                                                                                   
Out[12]: 
⎡          ⎛ ⎛ 2⎞    ⎞      ⎛ ⎛ 2⎞    ⎞           ⎛ 2⎞                   ⎛ ⎛ 2⎞    ⎞      ⎛ ⎛ 2⎞    ⎞           ⎛ 2⎞       ⎤
⎢          ⎜ ⎝t ⎠    ⎟      ⎜ ⎝t ⎠    ⎟           ⎝t ⎠                   ⎜ ⎝t ⎠    ⎟      ⎜ ⎝t ⎠    ⎟           ⎝t ⎠       ⎥
⎢          ⎜ℯ       1⎟      ⎜ℯ       1⎟   t   √π⋅ℯ    ⋅erf(t)            ⎜ℯ       1⎟      ⎜ℯ       1⎟   t   √π⋅ℯ    ⋅erf(t)⎥
⎢x(t) = C₁⋅⎜───── - ─⎟ + C₂⋅⎜───── + ─⎟ + ─ + ───────────────, y(t) = C₁⋅⎜───── + ─⎟ + C₂⋅⎜───── - ─⎟ - ─ + ───────────────⎥
⎣          ⎝  2     2⎠      ⎝  2     2⎠   2          4                   ⎝  2     2⎠      ⎝  2     2⎠   2          4       ⎦

[oscarbenjamin]

It's also worth noting that in situations where we get a different result depending on whether we collect first for t or the constants a change of variables in the constants can simplify the solution further and make it independent of collection order e.g.:

In [96]: simpsol(sol, [t], symbols('C1:4'))                                                                                                                   
Out[96]: 
⎡                                             ⎛ 2⎞                                                    ⎛ 2⎞       ⎤
⎢                                  ⎛ 2⎞       ⎝t ⎠                                         ⎛ 2⎞       ⎝t ⎠       ⎥
⎢         C₁   C₂   t   ⎛C₁   C₂⎞  ⎝t ⎠   √π⋅ℯ    ⋅erf(t)         C₁   C₂   t   ⎛C₁   C₂⎞  ⎝t ⎠   √π⋅ℯ    ⋅erf(t)⎥
⎢x(t) = - ── + ── + ─ + ⎜── + ──⎟⋅ℯ     + ───────────────, y(t) = ── - ── - ─ + ⎜── + ──⎟⋅ℯ     + ───────────────⎥
⎣         2    2    2   ⎝2    2 ⎠                4                2    2    2   ⎝2    2 ⎠                4       ⎦

In [97]: sol = [s.subs({C1:C3+C4, C2:C3-C4}) for s in sol]                                                                                                    

In [98]: simpsol(sol, [t], symbols('C1:4'))                                                                                                                   
Out[98]: 
⎡                               ⎛ 2⎞                                        ⎛ 2⎞       ⎤
⎢           ⎛ 2⎞                ⎝t ⎠                    ⎛ 2⎞                ⎝t ⎠       ⎥
⎢           ⎝t ⎠        t   √π⋅ℯ    ⋅erf(t)             ⎝t ⎠        t   √π⋅ℯ    ⋅erf(t)⎥
⎢x(t) = C₃⋅ℯ     - C₄ + ─ + ───────────────, y(t) = C₃⋅ℯ     + C₄ - ─ + ───────────────⎥
⎣                       2          4                                2          4       ⎦

In [99]: simpsol(sol, symbols('C1:3'), [t])                                                                                                                   
Out[99]: 
⎡                               ⎛ 2⎞                                        ⎛ 2⎞       ⎤
⎢           ⎛ 2⎞                ⎝t ⎠                    ⎛ 2⎞                ⎝t ⎠       ⎥
⎢           ⎝t ⎠        t   √π⋅ℯ    ⋅erf(t)             ⎝t ⎠        t   √π⋅ℯ    ⋅erf(t)⎥
⎢x(t) = C₃⋅ℯ     - C₄ + ─ + ───────────────, y(t) = C₃⋅ℯ     + C₄ - ─ + ───────────────⎥
⎣                       2          4                                2          4       ⎦

For single ODEs dsolve has constantsimp to simplify the solution by changing the constants. We should (eventually) do the same thing in dsolve_system.

[mijo2]

I will open the PR for simplification soon. But I wanted to point out that I made a mistake while designing the last type of systems. I didn't change the non-homogeneous term. I will rectify this along with an example in the upcoming PR.

[mijo2]

@oscarbenjamin When it comes to doit option, should I disable it with the simplification PR or should I just work on simplification for now?

[mijo2]

For single ODEs dsolve has constantsimp to simplify the solution by changing the constants. We should (eventually) do the same thing in dsolve_system.

Agreed. Or maybe we can extend the constantsimp to handle systems of ODEs. For now, we have to decide the preference for simplification, like constants first or terms of t first? Should I analyze which one will lead to a smaller number of operations on an average and make the decision?

[oscarbenjamin]

When it comes to doit option, should I disable it with the simplification PR

I'm not sure what you mean.

I think that some of the example systems in the tests will take far too long if dsolve_system attempts to evaluate the integrals so we need the option to avoid the integrals as well.

For now, we have to decide the preference for simplification, like constants first or terms of t first?

I think that it should be in terms of t first so that the coefficients of that are what would be simplified by constantsimp.

[mijo2]

I think that it should be in terms of t first so that the coefficients of that are what would be simplified by constantsimp.

Ok, I would work on this for now. I tried your method for one of the bigger test cases, and it took a lot of time to simplify it. So, there are definitely cases where it would be for the best for an end-user to not call the simplification method. I think we need to explicitly add about this in the documentation.

[mijo2]

I think that some of the example systems in the tests will take far too long if dsolve_system attempts to evaluate the integrals so we need the option to avoid the integrals as well.

I have added that and for now, I am keeping the default value of doit option as True.