Fix failing checkodesol examples

[oscarbenjamin]

We need to figure out the cases in the ODE tests where checkodesol doesn't work.

The tests in the ode module use checkodesol to check that the tested results are correct. So most tests there are of the form

eq = f(x).diff(x)-x
sol = Eq(f(x), x**2/2 + C1)
assert dsolve(eq, f(x)) == sol
assert checkodesol(eq, sol) == (True, 0)

The idea is that dsolve should give the solution of the ODE and checkodesol should verify that solution. So if everything is good then both those asserts will work. You can see examples here:

sympy/sympy/solvers/tests/test_ode.py

Line 3106 in c6726a1

def test_nth_algebraic():

However checkodesol doesn't always work. In #15957 checkodesol was improved to work in more cases and added to most tests in test_ode.py. There are a number of cases where it still didn't work which are now marked as FIXME in the tests e.g.:

sympy/sympy/solvers/tests/test_ode.py

Line 1778 in c6726a1

# FIXME: assert checkodesol(eq, sol) == (True, [0]) # Hangs...

I listed a summary of these cases in this comment:
#15957 (comment)

What is needed now is to go through each case and try to work out whether the tested solution is correct or not. If the solution is correct then the fact that checkodesol doesn't work is a bug in checkodesol. If the tested solution is incorrect then it is a bug in dsolve for returning an incorrect solution.

[Teut2711]

Hmm it's quite a lot to understand.But I would wanna work on it.

[oscarbenjamin]

Go ahead. If you find an example where you know if the solution is correct or not then report it here and we can decide what the fix should be.

One way to do this would be to see what another CAS (e.g. Wolfram Alpha) thinks the solution should be.

[sharma-kunal]

@oscarbenjamin I was looking up to the issue and while checking for the equation

eq2 = (2*x*f(x) + 1)/f(x) + (f(x) - x)/f(x)**2*f(x).diff(x)

having this solution or not

sol2 = Eq(f(x), exp(C1 - x**2 + LambertW(-x*exp(-C1 + x**2))))

dsolve returned the above solution but using Wolfram Alpha gave

Is it that dsolve is returning the wrong solution
Also above FIXME, it's commented #5080 is blocking testing for the solution

[oscarbenjamin]

Adding the solution above into the tests would be a reasonable PR.

[oscarbenjamin]

Okay so checkodesol doesn't like the solution returned by dsolve:

In [61]: eq2 = (2*x*f(x) + 1)/f(x) + (f(x) - x)/f(x)**2*f(x).diff(x)                                                              

In [62]: sol2 = Eq(f(x), exp(C1 - x**2 + LambertW(-x*exp(-C1 + x**2))))                                                           

In [63]: assert dsolve(eq2) == sol2                                                                                               

In [64]: eq2                                                                                                                      
Out[64]: 
            d                      
(-x + f(x))⋅──(f(x))               
            dx         2⋅x⋅f(x) + 1
──────────────────── + ────────────
        2                  f(x)    
       f (x)                       

In [65]: sol2                                                                                                                     
Out[65]: 
                          ⎛           2⎞
              2           ⎜    -C₁ + x ⎟
        C₁ - x  + LambertW⎝-x⋅ℯ        ⎠
f(x) = ℯ                                

In [66]: checkodesol(eq2, sol2)                                                                                                   
Out[66]: 
(False,
 2*x**3*exp(-C1 + x**2 - LambertW(-x*exp(-C1 + x**2))) + 2*x**2*LambertW(-x*exp(-C1 + x**2)) + x*exp(-C1 + x**2 - LambertW(-x*exp(-C1 + x**2))) + LambertW(-x*exp(-C1 + x**2)))

The answer from Wolfram is different and checkodesol thinks it is correct:

In [67]: sol_wolfram = -x / LambertW(-exp(x**2+C1)*x)                                                                             

In [68]: sol_wolfram                                                                                                              
Out[68]: 
         -x          
─────────────────────
        ⎛          2⎞
        ⎜    C₁ + x ⎟
LambertW⎝-x⋅ℯ       ⎠

In [69]: checkodesol(eq2, sol_wolfram)                                                                                            
Out[69]: (True, 0)

The answer from Wolfram is definitely not equivalent under any change of constant or anything so... this is a bug in dsolve. More specifically it is a bug in the 1st_exact solver:

In [72]: classify_ode(eq2)                                                                                                        
Out[72]: ('1st_exact', '1st_power_series', 'lie_group', '1st_exact_Integral')

I've tried the lie_group solver (dsolve(eq2, hint='lie_group') but it raises NotImplementedError.

We can see where the answer comes from a little bit with:

In [73]: dsolve(eq2, hint='1st_exact_Integral')                                                                                   
Out[73]: 
⎛                 ⌠                         ⎞│           
⎜⌠                ⎮ ⎛  ⌠                ⎞   ⎟│           
⎜⎮ ⎛      1⎞      ⎮ ⎜  ⎮ -1       1   x ⎟   ⎟│           
⎜⎮ ⎜2⋅x + ─⎟ dx + ⎮ ⎜- ⎮ ─── dx + ─ - ──⎟ dy⎟│       = C₁
⎜⎮ ⎝      y⎠      ⎮ ⎜  ⎮   2      y    2⎟   ⎟│           
⎜⌡                ⎮ ⎜  ⎮  y           y ⎟   ⎟│           
⎜                 ⎮ ⎝  ⌡                ⎠   ⎟│           
⎝                 ⌡                         ⎠│y=f(x)     

In [74]: _.doit()                                                                                                                 
Out[74]: 
 2    x                   
x  + ──── + log(f(x)) = C₁
     f(x)                 

In [75]: solve(_, f(x))                                                                                                           
Out[75]: 
⎡                   ⎛           2⎞⎤
⎢       2           ⎜    -C₁ + x ⎟⎥
⎢ C₁ - x  + LambertW⎝-x⋅ℯ        ⎠⎥
⎣ℯ                                ⎦

So either 1st_exact is buggy or integrate/solve are going wrong somewhere.

Can you open a new issue describing this (or just linking to this comment)?

[Teut2711]

If I can get help then I'm interested in fixing this.

[oscarbenjamin]

Adding the solution above to the tests would be a reasonable PR.

[oscarbenjamin]

For fixing it you'll need to to debug what's going on in 1st exact. Somewhere in there is probably a bug. If you can follow through what happens in a debugger you might be able to find it.

First compare the integrals above with what Wolfram is showing. Maybe Wolfram can give you a step-by-step solution which would help with debugging.

[Teut2711]

Ya, I understand that.I will look into the module but first I' ll have to something called "sync it" since changed from desolve might not have incorporated.

[sharma-kunal]

@XtremeGood I was also trying to solve the issue. So I think you can solve this bug and I may try to solve the next bug (where checkodesol hangs). I think that must be fine with you.

[sharma-kunal]

@oscarbenjamin Can you please help me with the equation,

eq = f(x).diff(x, 5) + 11*f(x).diff(x) - 2*f(x)

It's solution as solved by dsolve is

C5*exp(x*CRootOf(_x**5 + 11*_x - 2, 0)) + (C1*sin(x*im(CRootOf(_x**5 + 11*_x - 2, 1))) + C2*cos(x*im(CRootOf(_x**5 + 11*_x - 2, 1))))*exp(x*re(CRootOf(_x**5 + 11*_x - 2, 1))) + (C3*sin(x*im(CRootOf(_x**5 + 11*_x - 2, 3))) + C4*cos(x*im(CRootOf(_x**5 + 11*_x - 2, 3))))*exp(x*re(CRootOf(_x**5 + 11*_x - 2, 3)))

So it was hanging previously so what I did was I assumed the roots of the equation to be constant and now the output is coming quite fast, but the problem is the solution given by Wolfram Alpha is

wolf_sol = C3*exp(x*CRootOf(x**5 + 11*x - 2, 3)) + C4*exp(x*CRootOf(x**5 + 11*x - 2, 4)) + C5*exp(x*CRootOf(x**5 + 11*x - 2, 0)) + C1*exp(x*CRootOf(x**5 + 11*x - 2, 1)) + C2*exp(x*CRootOf(x**5 + 11*x - 2, 2))

and both (wolf_sol and sol) are giving FALSE with checkodesol. I can't figure out how to proceed to the issue.

[oscarbenjamin]

@sharma-kunal the two solutions are equivalent. In the wolf_sol solution some of the CRootOfs are complex:

In [5]: [s.evalf(5) for s in solve(x**5 + 11*x - 2)]                                                                              
Out[5]: [0.1818, -1.3312 - 1.2895⋅ⅈ, -1.3312 + 1.2895⋅ⅈ, 1.2403 - 1.2901⋅ⅈ, 1.2403 + 1.2901⋅ⅈ]

So with a complex root a we can rewrite exp(a*x) as exp(re(a)*x)*(cos(im(a)*x) + I*sin(im(a)*x)). If we do that and rearrange the constants we can rewrite wolf_sol as the result returned by dsolve.

It isn't sufficient to assume that the roots are constant: the solution only works if they are roots of that polynomial.

The problem is that the algebraic properties of rootofs aren't fully implemented. See #15753 where I proposed one way to fix this.

[sharma-kunal]

@oscarbenjamin I can't figure out what @asmeurer meant in this comment #15753 (comment). Can you explain it ?

[oscarbenjamin]

See also #18377