In one paper I find these set of equations:
$$ u_1 = b( \omega_1^2 + \omega_2^2 + \omega_3^2 + \omega_4^2)$$ $$ u_2 = b(\omega_1^2 + \omega_2^2 - \omega_3^2 - \omega_4^2)$$ $$ u_3 = b(\omega_1^2 - \omega_2^2 + \omega_3^2 - \omega_4^2)$$ $$ u_4 = b(\omega_1^2 - \omega_2^2 - \omega_3^2 + \omega_4^2) $$
Where we put numerical values in $\omega_i$ and $\omega_i>0$ (the constraint is based on physical meaning) and you can obtain $u_i$ solving these equations.
But my objective is the inverse, I would like to give values to $u_i$ and obtain $\omega_i$. The problem is that I don't know the range or the specific values I need in $u_i$ to obtain meaningful (physically plausible) values of $\omega_i$ ($\omega_i>0$).
How could I calculate these values using Mathematica? Is there a library of python that I can use? Up until now I have seen examples of $Y=AX$ and obtain $Y$ with constraints in $X$, not obtain $Y$ with constraints in $Y$.
---edit: possible solution---
I don't know why I thought it was a harder problem. At least in Mathematica software is straight forward. The steps I follow was convert to matrix form ($U=AW$), inverse matrix to get $W=A^{-1}U$ and in the equations also write $w_i>0$ conditions (wolfram alpha solutions) (In the link I used $x,y,z,t$ instead of $\omega_i$).
The only problem now is to know numerically for a range of $\omega_i$ what range I have in $u_i$. But as seen in the solution, is not quite easy, it depends on the relations between different $u_i$.
Define $$v_i=\frac {u_i}b \qquad \text{and}\qquad x_i=\omega_i^2$$ and you face four linear equations for four unknowns.