SymPy "solves" a differential equation it shouldn't solve
Asked Answered
T

1

5

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?

Thain answered 8/1, 2018 at 10:48 Comment(0)
L
7

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))
Lafreniere answered 8/1, 2018 at 15:8 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.