SymPy dsolve returns different results for mathematically equivalent differential equations

Question

Here is the content of my script:

from sympy import *
x = symbols('x')
init_printing(use_unicode=True)

f = symbols('f', cls=Function)
diffeq = Eq(x**2 * f(x).diff(x, x) + x * f(x).diff(x) - f(x) , 1/((1+x**2)**(3)) )
print dsolve(diffeq, f(x))

This program returns the following output:

Eq(f(x), (C1*x**2 + C1 + C2*x**4 + C2*x**2 - 15*x**4*atan(x) - 15*x**3 - 18*x**2*atan(x) - 13*x - 3*atan(x))/(16*x*(x**2 + 1)))

But when I define the variable diffeq like this:

diffeq = Eq(f(x).diff(x, x) + f(x).diff(x)/x - f(x)/x**(2) , 1 / ((1+x**2)**(3) * x**(2)) )

then I receive the output:

Traceback (most recent call last):
  File "/home/foo/odeSympyTrial01.py", line 12, in <module>
print dsolve(diffeq, f(x))
  File "/usr/lib/python2.7/dist-packages/sympy/solvers/ode.py", line 625, in dsolve
    x0=x0, n=n, **kwargs)
  File "/usr/lib/python2.7/dist-packages/sympy/solvers/deutils.py", line 235, in _desolve
    raise NotImplementedError(dummy + "solve" + ": Cannot solve " + str(eq))
NotImplementedError: solve: Cannot solve Derivative(f(x), x, x) + Derivative(f(x), x)/x - f(x)/x**2 - 1/(x**2*(x**2 + 1)**3)

And when I define the variable diffeq like this:

diffeq = Eq(f(x).diff(x, x) * x**(2) + f(x).diff(x) * x**(2) /x - f(x) * x**(2) /x**(2) , 1* x**(2)/((1+x**2)**(3) * x**(2)) )

then I receive the output:

Eq(f(x), (C1*x**2 + C1 + C2*x**4 + C2*x**2 - 15*x**4*atan(x) - 15*x**3 - 18*x**2*atan(x) - 13*x - 3*atan(x))/(16*x*(x**2 + 1)))

In every one of these cases, the differential equations diffeq are mathematically equal. Therefore in my opinion dsolve() should return the same output for each case. Somebody please help me to understand why dsolve() returns an error in the second case. How should the nonhomogeneous linear ordinary differential equation be expressed to ensure dsolve() does not return an error?

Answer 1

Short explanation: the logic of SymPy ODE module is often naive and sometimes incorrect.

As written originally, with

x**2 * f(x).diff(x, x) + x * f(x).diff(x) - f(x)

this matches the form of Cauchy–Euler equation (also known as Euler's equation): the power of x in each coefficient is the order of derivative. SymPy detects this structure and applies an appropriate method. But if you divide by x**2,

f(x).diff(x, x) + f(x).diff(x)/x - f(x)/x**(2) 

this is no longer the case: the second derivative does not have the power of x**2 so the match fails. A more careful check could detect the latent Cauchy-Euler structure here, but that's not implemented, as one can see by looking at the source.

You can check that this is indeed what's going on with

classify_ode(diffeq, f(x))

which will return 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters' in the first case but not the second.

---

While looking at the source, one can also seen an example of wrong logic.

        if coeff.is_Mul:
            if coeff.has(f(x)):
                return False
            return x**order in coeff.args

For example, x**2*sin(x) will pass this check with order=2, which means that SymPy will mistake x**2*sin(x)*f(x).diff(x, x) - f(x) = 0 for Euler's equation. And indeed,

dsolve(x**2*sin(x)*f(x).diff(x, x) - f(x), f(x))

"solves" the equation incorrectly. Do not trust ODE solutions from SymPy.