1

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?

EricVonB
  • 238
  • 3
  • 11
  • I don't see what your 3rd example is getting at: the output is the same as in the first case, nothing unusual. –  May 19 '18 at 16:22
  • The differential equation solver in SymPy uses pattern matching, which can be sensitive to the exact form of an expression. Also sometimes more than one solver matches (you can see all the matching patterns with `classify_ode` and specify which one to use with the `hint` flag to `dsolve`). Each solver should produce a correct answer though, even if it is in a different form (the values of the constants may not be the same for different solutions). If you find a wrong answer, please open an issue for it in the SymPy issue tracker. – asmeurer May 22 '18 at 03:13

1 Answers1

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.