How to narrow down "Invalid argument types for subtraction" in a system of equations?

[Andreas Schuldei]

I have a system of longish equations and when calling .solve() on it, I get this traceback:

Traceback (most recent call last):
  File "C:/Users/Andreas Schuldei/PycharmProjects/lissajous-achse/hgü-kabel-detektion-symbolic.py", line 109, in <module>
    result = solve(equations, (CM0, theta_i, theta_j, theta_k, i, a, B_earth))
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\solvers\solvers.py", line 1096, in solve
    solution = _solve_system(f, symbols, **flags)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\solvers\solvers.py", line 1730, in _solve_system
    i, d = _invert(g, *symbols)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\solvers\solvers.py", line 3118, in _invert
    rhs -= indep
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\core\numbers.py", line 2194, in __sub__
    return Rational.__sub__(self, other)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\core\decorators.py", line 89, in __sympifyit_wrapper
    return func(a, b)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\core\numbers.py", line 1725, in __sub__
    return Number.__sub__(self, other)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\core\decorators.py", line 89, in __sympifyit_wrapper
    return func(a, b)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\core\numbers.py", line 733, in __sub__
    return AtomicExpr.__sub__(self, other)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\core\decorators.py", line 251, in _func
    return func(self, other)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\core\decorators.py", line 126, in binary_op_wrapper
    return f(self)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\core\decorators.py", line 127, in binary_op_wrapper
    return func(self, other)
  File "C:\Users\Andreas Schuldei\PycharmProjects\lissajous-achse\venv\lib\site-packages\sympy\vector\basisdependent.py", line 352, in __rsub__
    raise TypeError("Invalid argument types for subtraction")
TypeError: Invalid argument types for subtraction

Process finished with exit code 1

I stared at my code long and hard, and fixed all instances of wrong types that i could think of. Also, in order to catch the error earlier, I tried insterting .simplify() calls in places that might benefit from them - but the process ran for hours, doing the simplyfy() without even reaching the solve() call.

So I am asking for tricks to investigate the types and them fitting together, earlier. You can find my code attached. (please critique it, I want to learn!).

[Kalevi Suominen]

I do not know solve well but I suspect that it cannot handle vector equations. There should probably be separate equations for the components.

Kalevi Suominen

[Andreas Schuldei]

how can I split up the i, j, and k components of the equations? does coeff() work on equations, too? or do i need to split up each vector into .coeff(i, j, k) separately?

[Andreas Schuldei]

Obviously you looked at the code. Thank you! So you didn't see any obvious type screwups, either? Do you have other tips or remarks?

jks...@gmail.com schrieb am Samstag, 30. Oktober 2021 um 14:23:41 UTC+2:

[Andreas Schuldei]

Thank you, that seems to have been the issue. I split up the vectors into their coefficients and entered those vector components as separate equations. Now the python process is humming along at 1.8Gbyte memory and low to medium CPU utilization. How long is realistic to wait for a positive result?

jks...@gmail.com schrieb am Samstag, 30. Oktober 2021 um 14:23:41 UTC+2:

[Alan Bromborsky]

Sympy will not solve vector equations and the vector module seems to be an after thought.  Convert the vector equations to systems of non-linear scalar equations and hold your breath for the solver to give and answer -

https://docs.sympy.org/latest/modules/solvers/solvers.html

Part of the problem is that sympy is not a rule based computer algebra system.  Everything is hardcoded

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/3fd7fbd3-9dc9-47d0-9ed8-e1f92adfcfbfn%40googlegroups.com.

[David Bailey]

On 30/10/2021 16:06, Alan Bromborsky wrote:
> ("Invalid argument types for subtraction"

Wouldn't it be an idea if errors like this included the types of the
relevant expressions and/or copied the expressions into location where
they could be retrieved?

David

[Oscar Benjamin]

On Sat, 30 Oct 2021 at 15:41, Andreas Schuldei <and...@schuldei.org> wrote:
>
>
> Thank you, that seems to have been the issue. I split up the vectors into their coefficients and entered those vector components as separate equations. Now the python process is humming along at 1.8Gbyte memory and low to medium CPU utilization. How long is realistic to wait for a positive result?

Symbolic algorithms often have extremely bad algorithmic complexity so
it's hard to say.

I took a look at the equations you are solving from the code you
showed and I don't quite understand what the unknowns are e.g. the
first 4 equations are:

In [24]: equations[0]
Out[24]: (c₀₁ - c₁₁)⋅(c₀₁ - m₀₁) + (c₀₂ - c₁₂)⋅(c₀₂ - m₀₂) + (c₀₃ -
c₁₃)⋅(c₀₃ - m₀₃) = 0

In [25]: equations[1]
Out[25]: (c₀₁ - c₁₁)⋅(c₁₁ - m₁₁) + (c₀₂ - c₁₂)⋅(c₁₂ - m₁₂) + (c₀₃ -
c₁₃)⋅(c₁₃ - m₁₃) = 0

In [26]: equations[2]
Out[26]: (c₀₁ - c₂₁)⋅(c₂₁ - m₂₁) + (c₀₂ - c₂₂)⋅(c₂₂ - m₂₂) + (c₀₃ -
c₂₃)⋅(c₂₃ - m₂₃) = 0

In [46]: equations[3]
Out[46]: (c₀₁ - c₃₁)⋅(c₃₁ - m₃₁) + (c₀₂ - c₃₂)⋅(c₃₂ - m₃₂) + (c₀₃ -
c₃₃)⋅(c₃₃ - m₃₃) = 0

The first unknown listed in the call to solve is CM0:

In [27]: CM0
Out[27]: (c₀₁ - m₀₁) i_Sys_sensors + (c₀₂ - m₀₂) j_Sys_sensors + (c₀₃
- m₀₃) k_Sys_sensors

I don't know what it means to solve for this vector quantity. Which of
the symbols here are the unknowns (the c0i or the m0i)?

If you wanted the m0i then these equations are linear but if you
wanted the c0i then these first equations are a quadratic polynomial
system. Either way the complexity of the solutions is going to grow
quicker than you might expect.

Your next equation 4 equations are monstrous. Here's the first of them
(I guess this is actually 3 scalar equations once you take
components):

In [48]: print(equations[4])
Eq((-B01 + B_x + 2*a*i*mu_0*((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 -
m03)*sin(theta_i)*cos(theta_j))*((c01 - m01)*cos(theta_j)*cos(theta_k)
+ (c02 - m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))*cos(theta_j)*cos(theta_k)/(pi*(((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 - m03)*sin(theta_j))**2)**2) +
(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i))*(a*i*mu_0/(pi*(((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 - m03)*sin(theta_j))**2)) -
2*a*i*mu_0*((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))**2/(pi*(((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))**2)**2)))*Sys_sensors.i + (-B02 + B_y +
2*a*i*mu_0*((c01 - m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 -
m03)*sin(theta_i)*cos(theta_j))*((c01 - m01)*cos(theta_j)*cos(theta_k)
+ (c02 - m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))*sin(theta_k)*cos(theta_j)/(pi*(((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 - m03)*sin(theta_j))**2)**2) +
(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k))*(a*i*mu_0/(pi*(((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 - m03)*sin(theta_j))**2)) -
2*a*i*mu_0*((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))**2/(pi*(((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))**2)**2)))*Sys_sensors.j + (-B03 + B_z -
2*a*i*mu_0*((c01 - m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 -
m03)*sin(theta_i)*cos(theta_j))*((c01 - m01)*cos(theta_j)*cos(theta_k)
+ (c02 - m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))*sin(theta_j)/(pi*(((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 - m03)*sin(theta_j))**2)**2) +
(a*i*mu_0/(pi*(((c01 - m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 - m03)*sin(theta_j))**2)) -
2*a*i*mu_0*((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))**2/(pi*(((c01 -
m01)*(sin(theta_i)*sin(theta_j)*cos(theta_k) -
sin(theta_k)*cos(theta_i)) + (c02 -
m02)*(sin(theta_i)*sin(theta_j)*sin(theta_k) +
cos(theta_i)*cos(theta_k)) + (c03 - m03)*sin(theta_i)*cos(theta_j))**2
+ ((c01 - m01)*cos(theta_j)*cos(theta_k) + (c02 -
m02)*sin(theta_k)*cos(theta_j) - (c03 -
m03)*sin(theta_j))**2)**2))*sin(theta_i)*cos(theta_j))*Sys_sensors.k,
0)

If you do get an answer from solve then it will be extremely
complicated. I doubt though that solve will be able to give an answer
because it tries to solve polynomial systems in radicals and for a
fully symbolic system of equations like this you will hit up against
the Abel-Ruffini limit.

Probably you need to try a different basic approach to your problem
from the outset like maybe it's better to use quaternions rather than
rotation angles or something or maybe you shouldn't try to fully solve
the system algebraically.

The question is what did you want the solutions for? Presumably you
wanted to use them for some other purpose in which case maybe there's
another way to achieve that other purpose. If you just want to look at
the solutions and see what they are then I think that if you do manage
to get expressions for the solution of this system (from any CAS) they
are going to be too complicated to just "look at".

Oscar

[Andreas Schuldei]

Thank you very much, Oscar. That is excellent feedback.

You are right. The solver didn't find a solution after about 18h of running. For the solver to work, I had broken up my vector equation into vector component equations, and while doing that, I chose slightly different variables to solve for, too. The equations were still monstrous, though, and it is little wonder that the solver found no solution. I will give your suggestion to use quaternions a try once I understand what they are. :-)

I don't care if the solution I get is too complicated to look at, as long as I can implement it and feed it with my known values. I need it to spit out the one result I want - that vector CM0 (or CM1, CM2, or CM3 - in real life, those will be roughly similar), but even the value c02, c12, c22, or c32 would be equally good to know. I still refuse to throw a numerical search algorithm at this problem, even though that most likely would work, because I hope to better understand the problem (and solution) by solving it symbolically.

[Oscar Benjamin]

On Sun, 31 Oct 2021 at 20:10, Andreas Schuldei <and...@schuldei.org> wrote:
>
> Thank you very much, Oscar. That is excellent feedback.
>
> You are right. The solver didn't find a solution after about 18h of running. For the solver to work, I had broken up my vector equation into vector component equations, and while doing that, I chose slightly different variables to solve for, too. The equations were still monstrous, though, and it is little wonder that the solver found no solution. I will give your suggestion to use quaternions a try once I understand what they are. :-)

They might not help but I have sometimes found that problems involving
fully arbitrary 3D rotations can be handled more simply with
quaternions.

> I don't care if the solution I get is too complicated to look at, as long as I can implement it and feed it with my known values. I need it to spit out the one result I want - that vector CM0 (or CM1, CM2, or CM3 - in real life, those will be roughly similar), but even the value c02, c12, c22, or c32 would be equally good to know.

It might be possible to solve a subset of your equations for a subset
of your unknowns to get a partial analytic solution.

> I still refuse to throw a numerical search algorithm at this problem, even though that most likely would work, because I hope to better understand the problem (and solution) by solving it symbolically.

If your ultimate goal is to feed in values then why is it better to do
that after calling solve rather than before? You say that you want to
better understand the problem and solution but I think that if you got
a solution in general symbols here then it would be so monstrously
complicated that it wouldn't help to understand anything.

Analytic solutions are nice when they are possible and not overly
complicated. Numerics are more widely used than symbolics precisely
because many problems do not have useful analytic solutions but can be
handled numerically.


Oscar

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/e36b540d-6b4c-4471-9f9c-c09bdddf3ecfn%40googlegroups.com.

[Alan Bromborsky]

Also note that for 3D quaternions are subalgebra of the geometric
algebra.  Hamilton did not come up with the pseudoscalar which allows
scalars, vectors, quaternions, and the pseudoscalar to form a complete
algebra.  In Hamiltons quaternions you can add scalars to the basis
quaternions but you cannot add vectors to either. Remember Hamilton's i,
j, and k are isomorphic to the basis vectors but they are not vectors
themselves.

[Andreas Schuldei]

Alan, I don't understand everything you write here due to my lack of math understanding. Can you possibly explain it like I was five?

Do I understand correctly that quaternions are not the solution to my unique needs regarding rotations because quaternions only rotate around an axis through the origin? And  I have two rotations, and at least one of them is away from the origin, which would require an additional translation, which is incompatible with quaternions?

[Andreas Schuldei]

The technical advantage of a symbolic solution is shorter, constant execution time compared to a numeric search. If all goes well, this algorithm will be used on embedded systems, perhaps under CPU constraints.
.

[Oscar Benjamin]

On Mon, 1 Nov 2021 at 17:17, Andreas Schuldei <and...@schuldei.org> wrote:
> .
> Oscar schrieb am Sonntag, 31. Oktober 2021 um 23:50:57 UTC+1:
>>
>> On Sun, 31 Oct 2021 at 20:10, Andreas Schuldei <and...@schuldei.org> wrote:
>> >
>> > Thank you very much, Oscar. That is excellent feedback.
>> >
>> > You are right. The solver didn't find a solution after about 18h of running. For the solver to work, I had broken up my vector equation into vector component equations, and while doing that, I chose slightly different variables to solve for, too. The equations were still monstrous, though, and it is little wonder that the solver found no solution. I will give your suggestion to use quaternions a try once I understand what they are. :-)
>>
>> They might not help but I have sometimes found that problems involving
>> fully arbitrary 3D rotations can be handled more simply with
>> quaternions.
>>
>> > I don't care if the solution I get is too complicated to look at, as long as I can implement it and feed it with my known values. I need it to spit out the one result I want - that vector CM0 (or CM1, CM2, or CM3 - in real life, those will be roughly similar), but even the value c02, c12, c22, or c32 would be equally good to know.
>>
>> It might be possible to solve a subset of your equations for a subset
>> of your unknowns to get a partial analytic solution.
>>
>> > I still refuse to throw a numerical search algorithm at this problem, even though that most likely would work, because I hope to better understand the problem (and solution) by solving it symbolically.
>>
>> If your ultimate goal is to feed in values then why is it better to do
>> that after calling solve rather than before? You say that you want to
>> better understand the problem and solution but I think that if you got
>> a solution in general symbols here then it would be so monstrously
>> complicated that it wouldn't help to understand anything.
>>
>> Analytic solutions are nice when they are possible and not overly
>> complicated. Numerics are more widely used than symbolics precisely
>> because many problems do not have useful analytic solutions but can be
>> handled numerically.
>
> The technical advantage of a symbolic solution is shorter, constant execution time compared to a numeric search. If all goes well, this algorithm will be used on embedded systems, perhaps under CPU constraints.

It is not always the case that using a symbolic solution gives shorter
execution time. In some cases the symbolic solution explodes in
complexity to the point that just substituting the values into the
solution is slower than solving the original system numerically.

Here's a demonstration:

In [1]: a, b, c, d, e, x = symbols('a:e, x')

In [2]: p = a*x**4 + b*x**3 + c*x**2 + d*x + e

In [3]: values = {a:1,b:-2,c:3,d:-4,e:-5}

In [4]: r1, r2, r3, r4 = roots(p, x)

In [5]: %time result = r1.evalf(subs=values) # using symbolic
solution for one root
CPU times: user 296 ms, sys: 0 ns, total: 296 ms
Wall time: 297 ms

In [6]: %time results = nroots(p.subs(values)) # solve numerically
for all roots
CPU times: user 32 ms, sys: 0 ns, total: 32 ms
Wall time: 34.4 ms

If you print out r1 then you can see what the expression looks like
(it's the quartic formula). The formula is so large that we can
actually solve for the roots numerically more efficiently than by
evaluating the formula.

Of course both approaches here are quite slow because they are
implemented in Python and are not that well optimised. Also evalf and
nroots are not a fair comparison to standard numerical routines
because they both use multiprecision arithmetic for highly accurate
results. We can try using lambdify to get a more accurate timing for
pure Python code to evaluate this result:

In [9]: f = lambdify((a,b,c,d,e), r1, 'numpy')

In [10]: %time result = f(1,-2,3,-4,-5)
<string>:2: RuntimeWarning: invalid value encountered in sqrt
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 823 µs

In [11]: result
Out[11]: array(nan+0.j)

Although this is a lot faster it demonstrates another problem: this
just returns nan because the expression for the root cannot be
computed without complex numbers even though the root itself is real.

By contrast the numpy numerical roots function gives reasonably
accurate results very quickly:

In [12]: %time results = np.roots([1,-2,3,-4,-5])
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 663 µs

I think for your problem I would probably use sympy to:

1. Find a partial solution for some of the simpler equations.
2. Derive and simplify/optimise expressions for the equations and
their Jacobian (or perhaps its inverse)
3. Use code generation to implement these things in the target language
4. Run a numerical solver on the embedded system

Another possibility is that you could use sympy to get an approximate
solution. Being able to initialise a numerical routine with a good
guess makes it more robust but also means that in many applications
only a single Newton iteration is needed to get the target accuracy
resulting in a constant time algorithm.

--
Oscar