model2.showODEs(method='vanKampen') -> IndexError: list index out of range

[willfurnass]

Encountered when using Python 3.6 or 3.7 and latest acceptable versions of dependencies: see CI runs for #376

[willfurnass]

Reproduced with sympy 1.5.x but not 1.4.y.

[willfurnass]

With sympy 1.4 these lines after line 2893 in models.py:

                    ODE_1 = solve(ode1, Derivative(PhiDict[nvec[0]], t), dict=True)
                    ODE_2 = solve(ode2, Derivative(PhiDict[nvec[1]], t), dict=True)

result in [{Derivative(Phi_A, t): -Phi_A**2*r_{A} - Phi_A*Phi_B*r_{A} - Phi_A*Phi_B*s - Phi_A*a_{A} - Phi_A*g_{A} + Phi_A*r_{A} - Phi_B*g_{A} + g_{A}}] and [{Derivative(Phi_B, t): -Phi_A*Phi_B*r_{B} - Phi_A*Phi_B*s - Phi_A*g_{B} - Phi_B**2*r_{B} - Phi_B*a_{B} - Phi_B*g_{B} + Phi_B*r_{B} + g_{B}}] being assigned when running model2.showODEs(method='vanKampen') in the user guide Notebook.

But with Sympy 1.5.1 only empty lists are returned from the solve() calls.

This appears to be due to

Vlist_lhs, Vlist_rhs, substring = _get_orderedLists_vKE(stoich)

assigning different values to Vlist_lhs for the two sympy versions.

(The latest docs for Sympy indicate that both the API and implementation of solve() are a bit of a mess and people should switch to using solveset() instead: https://docs.sympy.org/latest/modules/solvers/solveset.html)

[willfurnass]

Exploring further, with Sympy 1.4 simplify(-sqrt(\overline{V})*Derivative(Phi_A, t)*Derivative(P(eta_A, eta_B, eta_U, t), eta_A)) appears to just return its argument, whereas with SymPy 1.5.0 and 1.5.1 simplify returns 0 (sympy.core.numbers.Zero).

@jarmarshall Can you identify anything in the SymPy 1.5 or 1.5.1 release notes that might explain this?

[jarmarshall]

I don't see how that could ever be simplified to 0 with the assumptions SymPy has access to, so think it must be a bug...

The following (fairly) minimal working example recreates the problem for me in terms of different behaviour between sympy v1.4 and 1.5, so I will raise an issue on sympy's GitHub:

from sympy import simplify, sqrt, symbols, Symbol, Derivative, Function
V = Symbol(r'\overline{V}', real=True, constant=True)
ea = Symbol('eta_A')
eb = Symbol('eta_B')
eu = Symbol('eta_U')
pa = Symbol('Phi_A')
t = symbols('t')
P = Function('P')
simplify(expr = -sqrt(V)*Derivative(pa, t)*Derivative(P(ea, eb, eu, t), ea))

[jarmarshall]

Hmmm... checking a bit further, it's the simplification of the derivative of pa wrt t that returns 0 - presumably it treated the derivative as something it knew nothing about before, whereas now it assumes it is of the form pa*t^0 - need to check if this is a documented change before raising the issue...

[jarmarshall]

OK, my thinking is evolving on this - since the simplification is apparently happening inside a call to solve we need to figure out what we should actually be doing now - maybe that is an option for solve, maybe it is replacing with solveset, or maybe it is actually switching to dsolve...
I think this is our only remaining upper-bounded package dependency, and 1.5.* is pretty recent so 1.4 should be good for a while longer - are there any other downsides with sticking with 1.4 for a while longer until this is figured out?

[willfurnass]

Not that I know of. Given SymPy 1.5 was released in December 2019 I agree that pinning sympy <1.5 makes sense for the pending MuMoT release.

The task of comparing solve / dsolve / solveset and deciding on a way forward would definitely be suited to someone who's worked with SymPy much more than I have.