Add Riccati ODE solver to dsolve benchmarks

[naveensaigit]

See https://github.com/sympy/sympy/pull/21459

[naveensaigit]

@oscarbenjamin Some existing tests are failing. I fixed the failure in test_solve.py, but I'm unable to understand the dsolve systems test which is failing. Could you please take a look?

[oscarbenjamin]

Are you asking about this?

    def test_make_ode_01():
        ode, params = _make_ode_01()
        t, y, y0, k = params
        result = dsolve(ode, y[1](t))
        eq_assumption = sympy.Q.is_true(Ne(k[1], k[0]))
        refined = refine(result, eq_assumption)
        ignore = k + y0 + (t,)
>       int_const = [fs for fs in refined.free_symbols if fs not in ignore][0]
E       AttributeError: 'list' object has no attribute 'free_symbols'

benchmarks/tests/test_dsolve.py:15: AttributeError

I think that what has happened is that dsolve now returns a list where it previously returned a single Eq.

[naveensaigit]

Are you asking about this?

Yes. Now dsolve returns two solutions and both of them are Piecewise. The expected solution doesn't seem to be equal to either of them.

[oscarbenjamin]

The dsolve change is bisected to https://github.com/sympy/sympy/commit/6318773a65e9a9f5e1b06890a658e9e1b63eaecc from https://github.com/sympy/sympy/pull/21124

The previous output was:

In [1]: v = Function('v')

In [2]: t, k0, k1, u0 = symbols('t, k0, k1, u0')

In [3]: eq = Eq(Derivative(v(t), t), k0*u0*exp(-k0*t) - k1*v(t))

In [4]: eq
Out[4]: 
d                 -k₀⋅t          
──(v(t)) = k₀⋅u₀⋅ℯ      - k₁⋅v(t)
dt                               

In [5]: dsolve(eq)
Out[5]: 
       ⎛           ⎛⎧        k₁⋅t                    ⎞⎞       
       ⎜           ⎜⎪      -ℯ                        ⎟⎟       
       ⎜           ⎜⎪───────────────────  for k₀ ≠ k₁⎟⎟  -k₁⋅t
v(t) = ⎜C₁ + k₀⋅u₀⋅⎜⎨    k₀⋅t       k₀⋅t             ⎟⎟⋅ℯ     
       ⎜           ⎜⎪k₀⋅ℯ     - k₁⋅ℯ                 ⎟⎟       
       ⎜           ⎜⎪                                ⎟⎟       
       ⎝           ⎝⎩         t            otherwise ⎠⎠ 

The new output is

In [5]: dsolve(eq)
Out[5]: 
⎡                                                        ⎧       C₁⋅k₀                 C₁⋅k₁                 k₀⋅u₀                    ⎤
⎢       ⎧    -k₀⋅t            -k₀⋅t                      ⎪─────────────────── - ─────────────────── - ───────────────────  for k₀ ≠ k₁⎥
⎢       ⎪C₁⋅ℯ      + k₀⋅t⋅u₀⋅ℯ       for k₀ = k₁         ⎪    k₁⋅t       k₁⋅t       k₁⋅t       k₁⋅t       k₀⋅t       k₀⋅t             ⎥
⎢v(t) = ⎨                                       , v(t) = ⎨k₀⋅ℯ     - k₁⋅ℯ       k₀⋅ℯ     - k₁⋅ℯ       k₀⋅ℯ     - k₁⋅ℯ                 ⎥
⎢       ⎪           nan               otherwise          ⎪                                                                            ⎥
⎢       ⎩                                                ⎪                              nan                                 otherwise ⎥
⎣                                                        ⎩                                                                            ⎦

This is really just a convoluted way of writing the same thing. If k0 == k1 then the first solution is equal to the corresponding case from the Piecewise and the second solution is nan. If k0 != k1 then it's the other way round. They might look different but they are equivalent e.g.:

In [15]: sol1, sol2 = dsolve(eq)

In [16]: sol2.rhs.args[0][0]
Out[16]: 
       C₁⋅k₀                 C₁⋅k₁                 k₀⋅u₀       
─────────────────── - ─────────────────── - ───────────────────
    k₁⋅t       k₁⋅t       k₁⋅t       k₁⋅t       k₀⋅t       k₀⋅t
k₀⋅ℯ     - k₁⋅ℯ       k₀⋅ℯ     - k₁⋅ℯ       k₀⋅ℯ     - k₁⋅ℯ    

In [17]: sol2.rhs.args[0][0].collect(Symbol('C1'), cancel)
Out[17]: 
    -k₁⋅t          k₀⋅u₀       
C₁⋅ℯ      - ───────────────────
                k₀⋅t       k₀⋅t
            k₀⋅ℯ     - k₁⋅ℯ  

Of course the single Piecewise solution is the preferred form so this is a regression. Looks related to the comment here: https://github.com/sympy/sympy/pull/21124#discussion_r618720278

[oscarbenjamin]

I'm really not sure why that dsolve example is tested here in the benchmark test suite though.

[asmeurer]

I guess the benchmarks are tested to make sure they don't break, and also so you can run the tests against an older version of SymPy to ensure that the benchmarks actually do benchmark the right thing. The latter case is more important in my opinion. If we want to guard against regressions, we should add the tests to the SymPy test suite. And if the goal is just to be able to check validity of the benchmark, the tests here should be written in a much more generic way. Instead of asserting an exact solution, it should only check that it is mathematically correct, e.g., using checkodesol, or using some sort of numerical consistency checks for non-solver benchmarks where this isn't possible.