Search code examples
pythonsympysymbolic-mathodedifferential-equations

SymPy "solves" a differential equation it shouldn't solve


Here's what I did:

from sympy import *
x = symbols("x")
y = Function("y")
dsolve(diff(y(x),x) - y(x)**x)

The answer I get (SymPy 1.0) is:

Eq(y(x), (C1 - x*(x - 1))**(1/(-x + 1)))

But that's wrong. Both Mathematica and Maple can't solve this ODE. What's happening here?


Solution

  • A bug. SymPy thinks it's a Bernoulli equation

    y' = P(x) * y + Q(x) * y**n
    

    without checking that the exponent n is constant. So the solution is wrong.

    I raised an issue on SymPy tracker. It should be soon fixed in the development version of SymPy and subsequently in version 1.2. (As an aside, 1.0 is a bit old, many things have improved in 1.1.1 although not that one.)

    With the correction, SymPy recognizes there is no explicit solution and resorts to power series method, producing a few terms of the power series:

    Eq(y(x), x + x**2*log(C1)/2 + x**3*(log(C1)**2 + 2/C1)/6 + x**4*(log(C1)**3 + 9*log(C1)/C1 - 3/C1**2)/24 + x**5*(log(C1)**4 + 2*(log(C1) - 1/C1)*log(C1)/C1 + 2*(2*log(C1) - 1/C1)*log(C1)/C1 + 22*log(C1)**2/C1 - 20*log(C1)/C1**2 + 20/C1**2 + 8/C1**3)/120 + C1 + O(x**6))
    

    You don't have to wait for the patch to get this power series, it can be obtained by giving SymPy a "hint":

    dsolve(diff(y(x), x) - y(x)**x, hint='1st_power_series')
    

    Works better with an initial condition:

    dsolve(diff(y(x), x) - y(x)**x, ics={y(0): 1}, hint='1st_power_series')
    

    returns

    Eq(y(x), 1 + x + x**3/3 - x**4/8 + 7*x**5/30 + O(x**6))