model2.showODEs(method='vanKampen') -> IndexError: list index out of range
Opened this issue · 7 comments
Encountered when using Python 3.6 or 3.7 and latest acceptable versions of dependencies: see CI runs for #376
Reproduced with sympy 1.5.x but not 1.4.y.
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)
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?
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))
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...
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?
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.