Solving linear system with 6 equations and 22 unknowns for six of the unknowns

I am trying to find the solution for the M variables in the following system.

$$0 = C_{b} M^{b}_{x} - M^{a}_{x} k_{2a} + M^{a}_{y} \left(\omega - \omega_{a}\right)\\ 0 = C_{a} M^{a}_{x} - M^{b}_{x} k_{2b} + M^{b}_{y} \left(\omega - \omega_{b}\right)\\ 0 = C_{b} M^{b}_{y} + M^{a}_{x} \left(- \omega + \omega_{a}\right) - M^{a}_{y} k_{2a} - M^{a}_{z} \omega_{1}\\ 0 = C_{a} M^{a}_{y} + M^{b}_{x} \left(- \omega + \omega_{b}\right) - M^{b}_{y} k_{2b} - M^{b}_{z} \omega_{1}\\ 0 = C_{b} M^{b}_{z} + \frac{M^{a}_{0}}{T_{1a}} + M^{a}_{y} \omega_{1} - M^{a}_{z} k_{1a}\\ 0 = C_{a} M^{a}_{z} + \frac{M^{b}_{0}}{T_{1b}} + M^{b}_{y} \omega_{1} - M^{b}_{z} k_{1b}$$

i.e. I am solving for the 6 $M_x$, $M_y$, $M_z$ variables, but I am not solving for the $M^a_0$, $M^b_0$ variables.

I have tried using Sympy to solve this, but it has been running for over 17 hours, and has not yet returned a solution.

I am interested in a fast way to solve this.

Edit:

Click this link to have a look at my code.

Edit 2:

For anyone who's interested, @Doug Lipinski's solution was very fast. However, as he mentioned the solution was VERY large. Below is just one of the solutions:

$$M^{b}_{z} = \frac{\frac{M^{a}_{0} \left(C_{a} - \frac{\omega_{1}^{2} \left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{\left(- k_{2b} - \frac{\left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right) \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)} - \frac{\left(- \omega + \omega_{b}\right) \left(\omega - \omega_{b}\right)}{\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}}\right) \left(- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)\right)}\right)}{T_{1a} \left(\frac{\omega_{1}^{2} \left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right) \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{\left(- k_{2b} - \frac{\left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right) \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)} - \frac{\left(- \omega + \omega_{b}\right) \left(\omega - \omega_{b}\right)}{\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}}\right) \left(- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)\right)^{2}} + \frac{\omega_{1}^{2}}{- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)} - k_{1a}\right)} - \frac{M^{b}_{0}}{T_{1b}}}{\frac{\omega_{1}^{2}}{- k_{2b} - \frac{\left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right) \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)} - \frac{\left(- \omega + \omega_{b}\right) \left(\omega - \omega_{b}\right)}{\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}}} - k_{1b} - \frac{\left(C_{a} - \frac{\omega_{1}^{2} \left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{\left(- k_{2b} - \frac{\left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right) \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)} - \frac{\left(- \omega + \omega_{b}\right) \left(\omega - \omega_{b}\right)}{\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}}\right) \left(- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)\right)}\right) \left(C_{b} - \frac{\omega_{1}^{2} \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{\left(- k_{2b} - \frac{\left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right) \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)} - \frac{\left(- \omega + \omega_{b}\right) \left(\omega - \omega_{b}\right)}{\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}}\right) \left(- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)\right)}\right)}{\frac{\omega_{1}^{2} \left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right) \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{\left(- k_{2b} - \frac{\left(C_{a} - \frac{C_{a} \left(- \omega + \omega_{b}\right) \left(\omega - \omega_{a}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right) \left(C_{b} - \frac{C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{b}\right)}{k_{2a} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)}\right)}{- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)} - \frac{\left(- \omega + \omega_{b}\right) \left(\omega - \omega_{b}\right)}{\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}}\right) \left(- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)\right)^{2}} + \frac{\omega_{1}^{2}}{- \frac{C_{a} C_{b} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)}{k_{2a}^{2} \left(\frac{C_{a} C_{b}}{k_{2a}} - k_{2b}\right)} - k_{2a} + \frac{1}{k_{2a}} \left(- \omega + \omega_{a}\right) \left(\omega - \omega_{a}\right)} - k_{1a}}}$$

Note that the equation has been forced onto multiple lines by this site.

I am going to try using some of Sympy's simplify functions to see if it can be simplified at all.

• Can you post your code? This should be trivially solvable via symbolic Gaussian elimination since the equations are linear in the M variables. You could solve it by hand in no more than 15 minutes. – Doug Lipinski May 25 '14 at 18:17
• Also, there are 8 M variables $\{M_x^a,M_x^b,M_y^a,M_y^b,M_z^a,M_z^b,M_0^a,M_0^b\}$, which ones do you want to solve for? I guess you probably want to treat $\{M_0^a,M_0^b\}$ as parameters. – Doug Lipinski May 25 '14 at 18:22
• Sorry. You are correct. I don't want to solve for the $M^a_0$, $M^b_0$ parameters. – daviewales May 25 '14 at 23:17
• I have added a link to my code. – daviewales May 25 '14 at 23:24
• I've found an easy way to make this more manageable. Instead of solving symbolically, then substituting numbers into the symbolic solution, just substitute the numbers into the original equations, and then solve with the matrix. – daviewales May 28 '14 at 11:40

I suggest using the solve_linear_system_LU function in SymPy rather than just the more generic solve. This works fine for me, although the result is ridiculously long. In fact it's so long that trying to display it nicely via MathJax (used by if you init_printing()) in iPython locks up my machine. Obviously my comment about solving this by hand in 15 minutes was way off, that would take forever. Try this:

from sympy import *

omegaA, omegaB, omega1, omega = symbols('\omega_a \omega_b \omega_1 \omega')
Max, Mbx, May, Mby, Maz, Mbz = symbols('M^a_x M^b_x M^a_y M^b_y M^a_z M^b_z')
Ma0, Mb0 = symbols('M^a_0 M^b_0')
k1a, k1b, k2a, k2b = symbols('k_1a k_1b k_2a k_2b')
Ca, Cb = symbols('C_a C_b')
T1a, T1b = symbols('T_1a T_1b')

system = Matrix( ( ( -k2a, Cb, -(omegaA-omega), 0, 0, 0, 0 ),
( Ca, -k2b, 0, -(omegaB-omega), 0, 0, 0 ),
( (omegaA-omega), 0, -k2a, Cb, -omega1, 0, 0 ),
( 0, (omegaB-omega), Ca, -k2b, 0, -omega1, 0 ),
( 0, 0, omega1, 0, -k1a, Cb, -Ma0/T1a ),
( 0, 0, 0, omega1, Ca, -k1b, -Mb0/T1b ) ) )

solve_linear_system_LU(system, [Max, Mbx, May, Mby, Maz, Mbz])


I get an answer very quickly, but it's extremely long. You should also double check my code since I haven't gone through it too carefully.

If there's not a way to simplify the solution, it seems incredibly impractical to actually use. I would bet it's faster to solve the system than to evaluate the symbolic solution.

• Your matrix is correct, and you are solving for the correct variables, so I assume the solution is also correct. I wonder if my computer actually solved the equations quickly, then locked up on the MathJax for 25 hours... – daviewales May 26 '14 at 3:58
• There are a number of substitutions that I haven't shown here which will remove a couple of the variables above. For instance, k1a can be expressed in terms of T1a and Ca. – daviewales May 26 '14 at 3:59
• I guess this is matrix is in band structure. I think it would be suitable to use sparse representations and solvers here. – Tolga Birdal May 26 '14 at 8:22
• @tbirdal I doubt you'll get any worthwhile advantage from sparse solvers for a 6x6 system. Especially when solving symbolically. Solving the system is actually very fast. – Doug Lipinski May 26 '14 at 13:29
• Solving the system is very fast, but as I am discovering, expanding or simplifying even a single result is very slow. Fortunately I don't really need the equation to be pretty. I should be able to use it as is. Thanks, by the way @DougLipinski. – daviewales May 26 '14 at 13:58

Doug Lipisnki demonstrated the power of symbolic calculations (which I use quite often), but I would like to offer a purely numeric approach.

We have a system $A x = b$ where A:[6×22], b:[6×1] and x:[22×1]

We want to solve for 6 unknowns $x_U=(x1,x2,x3,..)$ so we split the vector $x$ into known's and unkown's using the orthogonal projection matrices $K$ and $U$ such that

$$x = K x_K + U x_U$$

The structure of $K$ is [22×16] with a single 1 on each column ,on the row representing the known $x$ value. $U$ is [22×6] with a single 1 on each column at the row representing the unkown $x$ value.

For example if $x=(x1,x2,x3)$ had 3 values only, with $x2$ unknown, then $$K = \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{bmatrix} \;\;\;U = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}$$

To solve for the unknowns we have

$$\left. A K x_K + A U x_U = b \right\} x_U = (A U)^{-1}(b - A K x_K)$$

In the end you get to solve the following linear system:

$$\begin{bmatrix} -k_{2a} & C_b & \omega-\omega_a & 0 & 0 & 0 \\ C_a &-k_{2b} & 0 & \omega-\omega_b & 0 & 0 \\ \omega_a-\omega & 0 & -k_{2a}& C_b &-\omega_1 \\ 0& \omega_b-\omega & C_a & -k_{2b} & 0 & -\omega_1 \\ 0 & 0 & \omega_1 & 0 & -k_{1a} & C_b \\ 0 & 0 & 0 & \omega_1 & C_a & -k_{1b} \end{bmatrix} \begin{pmatrix} M_x^a \\ M_x^b \\ M_y^a \\ M_y^b\\M_z^a \\M_z^b \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ - \frac{M_0^a}{T_{1a}} \\ - \frac{M_0^b}{T_{1b}} \end{pmatrix}$$

• The only problem I see here is that the original system is not linear in all of the parameters (e.g. $C_bM_x^b$ terms etc.) so this projection formulation will not work. You must choose variables and parameters from the start to have a linear system. If you do this, you end up with exactly the system I used above. If you mean that this system should then be solved numerically, I think that's a good option and is what I meant by "I would bet it's faster to solve the system [numerically] than to evaluate the symbolic solution." – Doug Lipinski May 28 '14 at 14:05
• Typically a non linear system of equations is handled with a Newton-Raphson scheme. The original posting shows linear equations though. – ja72 May 28 '14 at 14:20
• My point is that the original system is NOT linear in all 22 unknowns (there are terms like $C_bM_x^b$ which is not linear if both are unknowns). You said x:[22x1] implying x is the vector of all unknowns. In that case, the original system cannot be written as Ax=b where A and b contain only known constants, try it. To have a linear system, you must immediately choose to treat some of the unknowns as parameters so they can be included in A and b and then they should not be included in x. – Doug Lipinski May 28 '14 at 14:35
• Being a variable does not mean it is an unknown. The wording in the title is misleading. The way I see it, there are only 6 unknowns at any time, and the rest are variables whose value is known a priori to solving the system. – ja72 May 28 '14 at 14:37
• I agree about the title, but this is based on specifying parameters versus variables. Why is $x$ [22x1] in your answer? Either $x$ is [6x1] and the projections are entirely unnecessary or $x$ is [22x1] and the system is nonlinear. – Doug Lipinski May 28 '14 at 14:45