Solve 4 simultaneous linear equations

I'm trying to solve 4 simultaneous linear equations because the result will be used in some software I'm creating.

The equations take the following form:

$Aw + Bx + Cy - Bz = P$

$-Bw + Ax + By + Cz = 0$

$Cw - Bx + Dy + Bz = 0$

$Bw + Cx - By + Dz = 0$

The values of A, B, C, D and P can be considered constant, as they're determined by known program variables (which change at runtime), whereas w, x, y and z are determined by these equations.

Because it's for a program I'm writing, I can let the computer do some of the heavy lifting, so it's not necessary for me to express w,x,y and z all in terms of only constants. That is, I'm really trying to find

z = constant (expressed in A,B,C,D & P)

y = F(z)

x = F(y,z)

w = F(x,y,z)

I made an attempt and found a solution which I will provide below, but in the larger context of the program I'm writing it isn't behaving as I expected and my reason for posting here is to determine if that's because of a problem with my math or if the problem lies elsewhere in my code.

My solution:

(Some grouping first because this took pages and as explained above, I don't really need to fully solve it):

$E = A^2 - B^2$

$F = B^2 + C^2$

$G = PBF$

$H = B^2 F - ACF + CDE -B^2 E$

$I = BCE + BDE + BCF + ABF$

$J = (D + C)(AC + B^2 ) - (A + C)(C^2 + B^2)$

$K = (A + C)(BC + BD)$

Solving for z, y, x & w:

$z = \frac{GK }{IJ - HK}$

$y = \frac{PBF + (B^2 F - ACF + CDE - B^2 E) z}{BCE + BDE + BCF + ABF}$

$x = \frac{-Dz - Cz}{A + C}$

$w = \frac{Ax + By + Cz}{B}$

Any insights would be greatly appreciated because short of attempting to solve the whole thing instead on the off-chance made a mistake I'm not really sure what else to do.

• Why don't you just write in a matrix form and use the matrix inverse? – Matti P. May 22 '18 at 6:17
• The programming language in which I have to write this software can't do matrix math. – AnAllergyToAnalogy May 22 '18 at 6:23
• If you're simply solving $x, y, z$ and $w$, given the constants $A, B, C, D$ and $P$, your best shot is really the matrix method. You can use Gaussian elimination to determine the variables. If your programming language doesn't support it, you should program it there yourself. It's not a very long code. – Matti P. May 22 '18 at 6:29
• Which programming language is it? – Matti P. May 22 '18 at 6:30
• Actionscript 2.0 :( I'm looking into the Gaussian method now though... been a number of years since I did any matrix stuff – AnAllergyToAnalogy May 22 '18 at 6:33

In matrix form, the system of equations is $$\left[ \begin{array}{cccc} A & B & C & -B \\ -B & A & B & C \\ C & - B & D & B \\ B & C & -B & D \end{array} \right] \left\{ \begin{array}{c} w \\ x \\y \\ z \end{array} \right\} = \left[ \begin{array}{c} P \\ 0 \\ 0 \\ 0 \end{array} \right]$$ There is an algebraic solution to this equation, and it can be found with the matrix inverse: $$\left\{ \begin{array}{c} w \\ x \\y \\ z \end{array} \right\} = \left[ \begin{array}{cccc} A & B & C & -B \\ -B & A & B & C \\ C & - B & D & B \\ B & C & -B & D \end{array} \right]^{-1} \left[ \begin{array}{c} P \\ 0 \\ 0 \\ 0 \end{array} \right]$$ The inverse matrix is a bit tricky to compute, but according to Wolfram Alpha, the result is the matrix $\frac{1}{\text{det}}M$ so that the components of $M$ are

• $M_{0,0}= AB^2 + 2CB + DB^2 + AD^2 -C^2 D$
• $M_{0,1}= BC^2 + 2BDC + BD^2$
• $M_{0,2}= C^3 + 2B^2 C - ADC + AB^2 + B^2 D$
• $M_{0,3}= -BC^2 - ABC - BDC - ABD$

The other coefficients you need not worry about, because they are multiplied by zeros in the equation. The determinant is $$\text{det} = A^2 (B^2 +D^2 ) + AB^2 (2C + D) - 2AC^2 D +B^2 (2C + D)^2 + C^4$$ So in short, you can compute the unknowns by $$\left\{ \begin{array}{c} w \\ x \\ y\\z \end{array} \right\} = \frac{1}{\text{det}} \left\{ \begin{array}{c} M_{0,0} \\ M_{0,1} \\ M_{0,2} \\ M_{0,3} \end{array} \right\} P$$

• All the answer seem to consider the regular case (det different from 0 in the answer above) only. If that is not guaranteed, the Gauss algorithms seems to be unavoidable. – Jens Schwaiger May 22 '18 at 9:20

I did not check your calculations but I proceeded in a very similar manner. What I obtained is $$z=-\frac{B P (A+C) (C+D)}{A^2 \left(B^2+D^2\right)+2 A \left(B^2 (2 C+D)-C^2 D\right)+B^2 (2 C+D)^2+C^4}$$ $$w=-\frac{ A \left(B^2+D^2\right)+B^2 (2 C+D)-C^2 D}{B (A+C) (C+D)}\color{red}{z}$$ $$y=\frac{ \left(B^2+C^2\right)\color{red}{w}+B (C+D)}{B^2-C D}\color{red}{z}$$ $$x=-\frac{B (\color{red}{w}-\color{red}{y})+D \color{red}{z}}{C}$$

You can write your system $M u = b$ as $[M|b]$ and perform Gauss elimination to bring it to upper diagonal form and infer the solution set from this. You have unknowns $u = (w, x, y, z)$ and the augmented matrix form is: $$\left[ \begin{array}{rrrr|r} A & B & C & -B & P \\ -B & A & B & C & 0 \\ C & -B & D & B & 0 \\ B & C & -B & D & 0 \end{array} \right] \to \left[ \begin{array}{rrrr|r} A+C & 0 & C+D & 0 & P \\ B & C & -B & D & 0 \\ C & -B & D & B & 0 \\ 0 & A+C & 0 & C+D & 0 \\ \end{array} \right]$$ I see no further further simplicfication. It now already depends on the particular values of $A+C$, $B$ and $C$ what kind of solution (you will end up with none, one, or infinite many solutions) will be there.

For this you will have to go through the Gauss algorithm. The resulting matrix will show you how many linear independent equations there are and if they are consistent with the non-homogenous part in the right-most column.

By augmented RREF we have

$$\begin{vmatrix}A&B&C&-B&P\\-B&A&B&C&0\\C&-B&D&B&0\\B&C&-B&D&0 \end{vmatrix}\to\begin{vmatrix}A&B&C&-B&P\\-B&A&B&C&0\\A+C&0&C+D&0&P\\0&A+C&0&C+D&0 \end{vmatrix}$$

then we can simplify the derivation by

$$w=p+ky=\frac{P}{A+C}-\frac{C+D}{A+C}y,\quad x=kz=-\frac{C+D}{A+C}z$$

and we can solve just the first two equations for the unknowns $y$ and $z$.

• @mvw ops yes of course, I fix. Thanks – user May 22 '18 at 6:50