Success of linsolve depends on how the unknowns are sorted (sympy>=1.7)

[ulugris]

I am having a very strange issue with linsolve. After upgrading to a sympy version equal or higher than 1.7, the following code, which worked perfectly on sympy 1.6, gets stuck at solving the linear system:

from sympy import symbols, Matrix, sin, cos, solve, linsolve
from sympy.physics.mechanics import dynamicsymbols

r, l, M, m, I, Im, J = symbols('r l M m I I_m J')
L, R, g, Kt, Kv, b, t = symbols('L R g K_t K_v b t')

x, th, V, i, tm, wm = dynamicsymbols('x theta V i tau_m omega_m')

thd = th.diff(t)
ths = th.diff(t, 2)
xd = x.diff(t)
xs = x.diff(t, 2)

vb = Matrix([x + l*sin(th), l*cos(th)]).diff(t)

Tb = (M*vb.dot(vb) + J*thd**2)/2
Tw = (m*xd**2 + I*(xd/r)**2)/2
Vb = M*g*l*cos(th)
L = Tb + Tw - Vb

eqs = Matrix([L.diff(thd).diff(t) - L.diff(th) + tm,
              L.diff(xd).diff(t) - L.diff(x) - tm/r])

eqm = Kt*i - Im*wm.diff(t) - b*wm - tm
eqm = eqm.subs(i, (V - Kv*wm)/R)
eqm = eqm.subs(wm, xd/r - thd)
eqs = eqs.subs(tm, solve(eqm, tm)[0])

var = Matrix([th, x, thd, xd, ths, xs])
jac = eqs.jacobian(var)
e0 = eqs.subs({v: 0 for v in var})
j0 = jac.subs({v: 0 for v in var})
lsys = e0 + j0*var

# This gets stuck
thsxs = Matrix(linsolve(lsys, [ths, xs]).args[0])

However, if I reverse the order of the unknowns, like this:

# This works
xsths = Matrix(linsolve(lsys, [xs, ths]).args[0])

then everything works fine, regardless of the sympy version.

[oscarbenjamin]

The internal implementation of linsolve was switched to use DomainMatrix in #18844. That should make it faster for cases like this. However in this particular case the domain is found to be EX because construct_domain does not like mixing a function with its derivative:

In [17]: construct_domain([f(t)])[0]
Out[17]: ZZ[f(t)]

In [18]: construct_domain([f(t).diff(t)])[0]
Out[18]: ZZ[Derivative(f(t), t)]

In [19]: construct_domain([f(t).diff(t)+f(t)])[0]
Out[19]: EX

The EX domain can be very slow which is essentially why this becomes slow.

If we replace the functions and derivatives with separate symbols then everything works fine and we get a reasonably simple expression for the result:

In [20]: t0,t1,t2 = symbols('t0,t1,t2')

In [21]: x0,x1,x2 = symbols('x0,x1,x2')

In [22]: rep = {th:t0, th.diff(): t1, th.diff(t, 2):t2, x:x0, x.diff():x1, x.diff(t, 2):x2}

In [23]: lsys.subs(rep)
Out[23]: 
⎡Kₜ⋅V(t)                 ⎛Kₜ⋅Kᵥ    ⎞      ⎛            2⎞      ⎛  Kₜ⋅Kᵥ   b⎞      ⎛  Iₘ      ⎞⎤
⎢─────── - M⋅g⋅l⋅t₀ + t₁⋅⎜───── + b⎟ + t₂⋅⎝Iₘ + J + M⋅l ⎠ + x₁⋅⎜- ───── - ─⎟ + x₂⋅⎜- ── + M⋅l⎟⎥
⎢   R                    ⎝  R      ⎠                           ⎝   R⋅r    r⎠      ⎝  r       ⎠⎥
⎢                                                                                             ⎥
⎢                                                           ⎛Kₜ⋅Kᵥ    ⎞      ⎛  Kₜ⋅Kᵥ   b⎞    ⎥
⎢                                                        t₁⋅⎜───── + b⎟   x₁⋅⎜- ───── - ─⎟    ⎥
⎢     Kₜ⋅V(t)      ⎛  Iₘ      ⎞      ⎛I    Iₘ        ⎞      ⎝  R      ⎠      ⎝   R⋅r    r⎠    ⎥
⎢   - ─────── + t₂⋅⎜- ── + M⋅l⎟ + x₂⋅⎜── + ── + M + m⎟ - ────────────── - ────────────────    ⎥
⎢       R⋅r        ⎝  r       ⎠      ⎜ 2    2        ⎟         r                 r            ⎥
⎣                                    ⎝r    r         ⎠                                        ⎦

In [24]: %time ok = linsolve(lsys.subs(rep), [t2, x2])
CPU times: user 120 ms, sys: 0 ns, total: 120 ms
Wall time: 118 ms

In [25]: %time ok = linsolve(lsys.subs(rep), [x2, t2])
CPU times: user 112 ms, sys: 0 ns, total: 112 ms
Wall time: 115 ms

I'm not sure whether it should be construct_domain or linsolve that should perform this conversion.

When you say that this worked perfectly with sympy 1.6 what exactly do you mean? When I try this with 1.6 I just get the empty set:

In [3]: linsolve(lsys, [xs, ths])
Out[3]: ∅

[ulugris]

Thank you for your quick answer, @oscarbenjamin!

linsolve(lsys, [xs, ths]) yields the correct solution in 1.6.2 and 1.6.4, I didn't actually test 1.6 itself. Maybe that's why you're getting an empty set.

One of the things I like about sympy over Matlab is being able to work with derivatives. I was used to creating intermediate symbols in Matlab, but I really like not needing to do so in sympy. Maybe this is not good practice anyway, and I should stick to using intermediate symbols.

[oscarbenjamin]

This is fixable. I'm just not sure whether it should be construct_domain that treats the derivatives as independent generators and constructs a ring like ZZ[x(t), Derivative(x(t), t)]. If that was changed then it would effect all Poly operations which is maybe a good thing although it might have unintended consequences somewhere.

Alternatively there could be a more targeted fix just for linsolve that would substitute particular things for symbols so that the polynomial code only needs to work with well-defined symbols as generators.

Possibly it does not make sense to do this in construct_domain because in general x(t) and Derivative(x(t), t) are not algebraically independent. For example if x(t) = sin(t) then x'(t) = cos(t) and we have that x(t)**2 + x'(t)**2 = 1. So even if we can treat x(t) as being "transcendental" I'm not sure we can treat that x(t) and x'(t) as simply being independent generators in general. If linsolve wants to do that then we can make something in linsolve that masks off all relevant functions with dummy symbols.

I'm honestly not sure though. @jksuom what do you think?

[oscarbenjamin]

But then I suppose x(t) could equally just be equal to zero or something so in so far as it makes sense for Poly to treat is as transcendental maybe it's also fine to treat x(t) and x'(t) as being independent generators...