dsolve doesn't know order-reducing substitutions

[oscarbenjamin]

If an ODE only involves derivatives of the dependent variable then it is possible to use a substitution to to reduce the order of the ODE. This usually makes it easier to solved but dsolve doesn't know how to do it yet.

Here's an example:

In [1]: x = Symbol('x')                                                                                                                        

In [2]: f = Function('f')                                                                                                                      

In [3]: eq = f(x).diff(x, 2) + x * f(x).diff(x)**2                                                                                             

In [4]: eq                                                                                                                                     
Out[4]: 
            2     2      
  ⎛d       ⎞     d       
x⋅⎜──(f(x))⎟  + ───(f(x))
  ⎝dx      ⎠      2      
                dx       

In [5]: dsolve(eq)                                                                                                                             
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)

Since f(x) only appears as a derivative we can use the substitution g(x) = f(x).diff(x) to obtain a first order ODE:

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

In [7]: eq_g = g(x).diff(x, 1) + x * g(x)**2                                                                                                   

In [8]: eq_g                                                                                                                                   
Out[8]: 
   2      d       
x⋅g (x) + ──(g(x))
          dx      

In [9]: dsolve(eq_g)                                                                                                                           
Out[9]: 
          2   
g(x) = ───────
             2
       C₁ + x 

This is a solution of the original ODE since

In [10]: sol = dsolve(eq_g).rhs                                                                                                                

In [12]: solf = integrate(sol, x)                                                                                                              

In [13]: solf                                                                                                                                  
Out[13]: 
      _____    ⎛         _____    ⎞       _____    ⎛       _____    ⎞
     ╱ -1      ⎜        ╱ -1      ⎟      ╱ -1      ⎜      ╱ -1      ⎟
-   ╱  ─── ⋅log⎜- C₁⋅  ╱  ───  + x⎟ +   ╱  ─── ⋅log⎜C₁⋅  ╱  ───  + x⎟
  ╲╱    C₁     ⎝     ╲╱    C₁     ⎠   ╲╱    C₁     ⎝   ╲╱    C₁     ⎠

In [14]: simplify(solf.diff(x, 2) + x*solf.diff(x)**2)                                                                                         
Out[14]: 0

It would be good to add a new solver in sympy/solvers.ode.py. The idea would be that the solver matches ODEs that only include on derivatives of the dependent variable. It would generate the ODE by substitution

In [15]: eq                                                                                                                                    
Out[15]: 
            2     2      
  ⎛d       ⎞     d       
x⋅⎜──(f(x))⎟  + ───(f(x))
  ⎝dx      ⎠      2      
                dx       

In [16]: eq.subs(f(x).diff(x), g(x))                                                                                                           
Out[16]: 
   2      d       
x⋅g (x) + ──(g(x))
          dx 

Then the solver should recursively call dsolve to solve that ODE and (if that works) integrate the result to get the solution to the original ODE.

This solver should appear in the matching order after all of the methods for linear ODEs but before the nonlinear methods (e.g. Liouville) here :https://github.com/sympy/sympy/blob/master/sympy/solvers/ode.py#L309.

[Teut2711]

Hi, I 'm interested in fixing this bug but I' will have to see how to convert my math to symbolic computations.

[Teut2711]

Even if you try this eq using dsolve => f(x).diff(x, 1) + x * f(x).diff(x)**2=0 it gives only the constant solution but not the other one. Whereas solve gives the other solution that too in f'(x) form not f(x). Definitely needs a lot of fixes.

[oscarbenjamin]

Hi @XtremeGood. If you want to work on this you'll need to use the latest master version of SymPy. The equation you show has two solutions under current master:

In [2]: f(x).diff(x, 1) + x * f(x).diff(x)**2                                                                                                  
Out[2]: 
            2           
  ⎛d       ⎞    d       
x⋅⎜──(f(x))⎟  + ──(f(x))
  ⎝dx      ⎠    dx      

In [3]: dsolve(_)                                                                                                                              
Out[3]: [f(x) = C₁, f(x) = C₁ - log(x)]

In terms of how to code this the method could be roughly as simple as

eq = f(x).diff(x, 2) + x * f(x).diff(x)**2                                                                                             
new_eq = eq.subs(f(x).diff(x), g(x))                                                                                                   
if new_eq.has(f): # if new_eq has f then f appears outside of derivatives
    raise NotImplementedError
sol = dsolve(new_eq, g(x)) 
sol_f = Eq(f(x), integrate(sol.rhs, x))                                                                                               

The more complicated part is working out how to slot that into dsolve's machinery. It needs to be split up into a matcher and a solver. The idea is that the matcher looks at the ODE and says whether this method will work. The solver later actually makes the solution. As a starting point here the matcher can just check new_eq.has(f). That would need to go somewhere here:
https://github.com/sympy/sympy/blob/master/sympy/solvers/ode.py#L1362

Then the solver part has to be a function called ode_<something>, perhaps ode_nth_order_substitution.

[Teut2711]

Hi @oscarbenjamin
I will check that today.Thanks for the help. :)
The problem is that classify_ode() returns empty ( ),i.e. zero number of hints.Gotta do something about it
and that equation is a special case of Louville equation.
https://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Liouville

[Teut2711]

Hi @oscarbenjamin
I would add a method in dsolve which which recursively tries solving it for f(x). If not solvable then try substituting f'(x) and solves for f'(x) and then for f(x) upto order of equation according to what you are saying.Also I thing a condition should be set to try for splitting it into factors.
I still wonder why sympy is not solving the equation by "Liouville ode" tactic which is already implemented.

[oscarbenjamin]

Hi again @XtremeGood.

If you add the code to match this type of ODE where I suggested then classify_ode will not return empty.

Perhaps you should open a new issue about the Liouville solver. The type of ODE referred to by this issue is more general than Liouville: it could be of any order.

[Teut2711]

Hi @oscarbenjamin
I googled about the liouville equation but couldn't find any major content.It says it has some things to do with lie groups which is not even taught in my college till even 3rd year of math course .Nevertheless for the time being I will write the recursion till maybe tomorrow and the thing you are telling is quite necessary in solving many equations,About liouville equations, I will try studying Lie groups first and try to make it work.

[Teut2711]

Hi @oscarbenjamin
One thing caught my mind that if the equation contains f(x),f'(x),f''(x)etc then only the lowest order should be substituted .To make it more clear consider eq = f(x).diff(x, 2) + x * f(x).diff(x)**2 + f(x).Now if I substitute f'(x)=g(x) it would change g(x)=f'(x) but f(x) would remain the same.Now what should I think to do?Again it should give NotImplemented.

Check my github for the current code which would work for the function which we were expecting:
https://github.com/XtremeGood/Sympy-Contrib

[oscarbenjamin]

This method should not be used if f(x) appears outside of derivatives. In dsolve it's not actually necessary to return NotImplementedError. We just need to test for this when deciding if the method matches the ODE. I recommend just to test if eq.subs(...).has(f(x)).

[Teut2711]

Done it, please see,I do have a github account but 'm not familiar to merge,pull requests etc.
https://github.com/XtremeGood/Sympy-Contrib/blob/master/Sympy%20Edit.py

If you see it good can you please tell me how to ask them for changes.

[oscarbenjamin]

Okay that looks like a good start. What we need to do now is get it into the SymPy codebase and get it up as a pull request. For this you'll need to install git. Let's not worry right now about making it perfect. Once you have a pull request we can make comments on how to improve it and we can see what happens when running the tests. I'll describe the steps but you'll probably want to read some guides. Here's the SymPy guide to contributing: https://github.com/sympy/sympy/wiki/introduction-to-contributing

In terms of git the steps you need to follow are:

1. Fork the SymPy repository in github
2. Install git
3. Clone your fork with git clone
4. Create a branch for this work git checkout -b dsolve_order
5. Make some changes to the code (just add what you have now as a new function in ode.py, we can improve it later)
6. Commit those changes git commit
7. Push them to your fork git push -u origin dsolve_order
8. Open a pull request in github (make sure you link back to this issue #15881.

Then we can talk about the rewriting the code so it fits into dsolve.

[Teut2711]

Thanks a lot! I will try to do it as fast as I can.

[Teut2711]

Done it.I have opened a pull request named "Sympyedit1". Please see it.It is showing error in spacing.

[ShubhamKJha]

I would like to work on this issue. Is there any lead on the way it should be worked upon?

[oscarbenjamin]

This is mostly done in #15992

[Teut2711]

Am I a part of this merge?

[asmeurer]

Yes, if you look the first commit on that pull request is your commit.

[Teut2711]

Thanks @asmeurer, @oscarbenjamin and @smichr for the great help in first merge!!

[oscarbenjamin]

Thanks @XtremeGood for your help!

There are many more issue with dsolve if you look...