dsolve fails to integrate simple ODE

109 views
Skip to first unread message

Oscar Benjamin

unread,
Sep 11, 2018, 6:53:26 PM9/11/18
to sympy
Hi,

I'm not sure if I'm missing a trick but I've been trying to use dsolve and it seems it doesn't work in many simple cases. I've put some examples below, tested with sympy 1.2 installed using pip. I don't know if any of the below is me using dsolve incorrectly or should be considered a bug or is just a reflection of it being a work in progress. 

I want to solve the heat/diffusion equation for heat in a cylinder with an exponential heat term:

$ isympy

In [1]: eqn = Derivative(Derivative(f(x),x)*x,x)/x - exp(x)


In [2]: eqn

Out[2]: 

       d ⎛  d      

       ──⎜x⋅──(f(x))⎟

   x   dx⎝  dx     

- ℯ  + ──────────────

             x       


In [3]: dsolve(eqn, f(x))

---------------------------------------------------------------------------

NotImplementedError


It strikes me as odd that sympy can't solve this since it's essentially an algebraic rearrangement to get f here and all that is needed is to understand that d/dx can be integrated. Sympy can do the integrals:

In [16]: res = expand((simplify(eqn * x).integrate(x) + C1)/x).integrate(x) + C2

    ...: 


In [17]: res

Out[17]: 

                         x        

C₁⋅log(x) + C₂ + f(x) - ℯ  + Ei(x)


In [18]: solve(res, f(x))

Out[18]: 

                  x       

⎣-C₁⋅log(x) - C₂ + ℯ  - Ei(x)⎦


The algorithm for this is essentially the same as solving an algebraic equation. If I replace the derivatives with logs then solve can do it:

In [23]: eqn = log(log(f(x))*x)/x - exp(x)


In [24]: eqn

Out[24]: 

   x   log(x⋅log(f(x)))

- ℯ  + ────────────────

              x        


In [25]: solve(eqn, f(x))

Out[25]: 

    x⎤

  x⋅ℯ ⎥

⎢ ℯ   

⎢ ─────⎥

  x 

⎣ℯ    


If I use f(x).diff(x) rather than Derivative then it looks more complicated but should work:

In [26]: eqn = ((f(x)).diff(x)*x).diff(x)/x - exp(x)


In [27]: eqn

Out[27]: 

           2                 

          d          d       

       x⋅───(f(x)) + ──(f(x))

           2         dx      

   x     dx                  

- ℯ  + ──────────────────────

                 x           


In [28]: dsolve(eqn)

---------------------------------------------------------------------------

NotImplementedError


The automatic expansion of the product rule here leaves us in a slightly trickier position but again there is a general rule that sympy is overlooking here: Use a substitution f' = g to bring it to 1st order:

In [32]: eqn = (g(x)*x).diff(x)/x - exp(x)


In [33]: eqn

Out[33]: 

         d              

       x⋅──(g(x)) + g(x)

   x     dx             

- ℯ  + ─────────────────

               x        


In [34]: dsolve(eqn, g(x))

Out[34]: 

                  x

       C₁    x    

g(x) = ── + ℯ  - ──

       x         x 


Then strangely I get:

In [41]: dsolve(eqn, g(x)).integrate(x)

Out[41]: 

                             

            ⎮ ⎛           x⎞   

          ⎮ ⎜C₁    x   ℯ ⎟   

⎮ g(x) dx = ⎮ ⎜── + ℯ  - ──⎟ dx

          ⎮ ⎝x         x ⎠   

               


But if I just integrate the rhs it does the integral:

In [40]: dsolve(eqn, g(x)).rhs.integrate(x)

Out[40]: 

             x        

C₁⋅log(x) + ℯ  - Ei(x)


--
Oscar

Aaron Meurer

unread,
Sep 11, 2018, 7:06:54 PM9/11/18
to sy...@googlegroups.com
The algorithms to solve it aren't implemented. In factored form, the
equation can be solved by integrating as you mentioned. There is an
issue to implement this algorithm, but it hasn't been done yet
https://github.com/sympy/sympy/issues/6259.

In expanded form it is a Euler equation, but the algorithm isn't smart
enough to recognize it (you need to multiply through by x**2). If you
do that you get an answer:

>>> dsolve((eqn.doit()*x**2).expand())
Eq(f(x), C1 + C2*log(x) + (x - 1)*exp(x)*log(x) - Integral(x*exp(x)*log(x), x))

I opened https://github.com/sympy/sympy/issues/15217 for this.

I guess based on your other calculation the answer should be
expressible via Ei, which SymPy doesn't know how to do for that
integral.

Regarding the final issue, I'm not sure why integrate() on an Eq
doesn't evaluate the integral, but if you call doit() it computes it.

>>> dsolve(eqn, g(x)).integrate(x).doit()
Eq(Integral(g(x), x), C1*log(x) + exp(x) - Ei(x))

I opened https://github.com/sympy/sympy/issues/15218 for this.

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/CAHVvXxROOENnFvfstx2ODdTC51JTYfaXgZMz%3DinO8HDKP2pSTA%40mail.gmail.com.
> For more options, visit https://groups.google.com/d/optout.

Oscar Benjamin

unread,
Sep 11, 2018, 7:28:55 PM9/11/18
to sympy
Thanks for the quick response. I have one more while I'm here:

In [43]: eqn = exp(f(x).diff(x)-f(x))


In [44]: eqn

Out[44]: 

         d       

 -f(x) + ──(f(x))

         dx      

               


In [45]: dsolve(eqn, f(x))

---------------------------------------------------------------------------

IndexError                                Traceback (most recent call last)

<ipython-input-45-b99728060ab1> in <module>()

----> 1 dsolve(eqn, f(x))


~/current/sympy/venv/lib/python3.6/site-packages/sympy/solvers/ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)

    662             # The key 'hint' stores the hint needed to be solved for.

    663             hint = hints['hint']

--> 664             return _helper_simplify(eq, hint, hints, simplify, ics=ics)

    665 

    666 def _helper_simplify(eq, hint, match, simplify=True, ics=None, **kwargs):


~/current/sympy/venv/lib/python3.6/site-packages/sympy/solvers/ode.py in _helper_simplify(eq, hint, match, simplify, ics, **kwargs)

    687         # attempt to solve for func, and apply any other hint specific

    688         # simplifications

--> 689         sols = solvefunc(eq, func, order, match)

    690         if isinstance(sols, Expr):

    691             rv =  odesimp(sols, func, order, cons(sols), hint)


~/current/sympy/venv/lib/python3.6/site-packages/sympy/solvers/ode.py in ode_lie_group(eq, func, order, match)

   5460         else:

   5461             y = Dummy("y")

-> 5462             h = sol[0].subs(func, y)

   5463 

   5464     if xis is not None and etas is not None:


IndexError: list index out of range


Not sure what the right response is here but perhaps not an IndexError.

--
Oscar


Aaron Meurer

unread,
Sep 11, 2018, 7:30:32 PM9/11/18
to sy...@googlegroups.com
IndexError indicates a bug. Can you open an issue for it?

Aaron Meurer

On Tue, Sep 11, 2018 at 5:28 PM, Oscar Benjamin
> https://groups.google.com/d/msgid/sympy/CAHVvXxROL%2BquRBfLhRBP_DSzXpV7KE64A7KeQ%3D5VpXdrMnsCkQ%40mail.gmail.com.

Oscar Benjamin

unread,
Sep 11, 2018, 7:35:37 PM9/11/18
to sympy
Reply all
Reply to author
Forward
0 new messages