Solving System of linear equation using LU decomposition

linear algebralu decompositionmatricesmatrix decompositionsystems of equations

I am working on a simultaneous linear equation problem using LU decomposition
and I'm unsure if this is the correct approach/answer to solve a system of simultaneous equations using LU decomposition. I'm looking for a way to check my resulting calculations and to understand if I have done the calculations correctly. This is my first time doing LU decomposition or doing matrix calculations at this level.

$5×1 + 6×2 + 2.3×3 + 6×4 = 4$

$9×1 + 2×2 + 3.5×3 + 7×4 = 5$

$3.5×1 + 6×2 + 2×3 + 3×4 = 6.7$

$1.5×1 + 2×2 + 1.5×3 + 6×4 = 7.8$

\begin{bmatrix}5&6&2.3&6\\9&2&3.5&7\\3.5&6&2&3\\1.5&2&1.5&6\end{bmatrix}

I multiply the top, third and fourth row by 10 and the second by 2 to make it easier to work with.

\begin{bmatrix}50&60&23&60\\18&4&7&14\\35&60&20&30\\15&20&15&60\end{bmatrix}

I calculated the L matrix as:

\begin{bmatrix}1&0&0&0\\9/25&1&0&0\\7/10&45/44&1&0\\3/10&5/44&1.584572&1\end{bmatrix}

and U matrix as:
\begin{bmatrix}50&60&23&60\\0&88/25&-1.28&-38/5\\0&0&5.2&-93/22\\0&0&0&12.6\end{bmatrix}

I have y as the following:
\begin{bmatrix}4\\5\\6.7\\7.8\end{bmatrix}

and solved the system with the following values for x
\begin{bmatrix}4\\3.56\\0.2590\\5.734906\end{bmatrix}

Best Answer

Your answer has a few issues, it was close.

You made some sign errors in values of $L$, and that messed up the final $L$.

For example, you should have $−45/44$ and $−5/44$ and that looks like it caused $L_{43}$ to be wrong. Of course, these issues caused similar errors in $U$.

The result should be

$$A = \begin{pmatrix}50&60&23&60\\18&4&7&14\\35&60&20&30\\15&20&15&60\end{pmatrix} = LU = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ \frac{9}{25} & 1 & 0 & 0 \\ \frac{7}{10} & -\frac{45}{44} & 1 & 0 \\ \frac{3}{10} & -\frac{5}{44} & \frac{175}{57} & 1 \\ \end{array} \right)\left( \begin{array}{cccc} 50 & 60 & 23 & 60 \\ 0 & -\frac{88}{5} & -\frac{32}{25} & -\frac{38}{5} \\ 0 & 0 & \frac{57}{22} & -\frac{435}{22} \\ 0 & 0 & 0 & \frac{1935}{19} \\ \end{array} \right)$$

Hopefully, you can take it from here.

Notes:

$1.$ You can always verify your result by substituting the values you got for $x$ back into the original system and seeing if the LHS = RHS.
$2.$ You could have also multiplied $LU$ and see if that returned $A$.
$3.$ Also, please recall that you multiplied your equations by $10$ and $2$, but forgot to make the corresponding update to $y$ for the final steps.

We now want to use back substitution to solve

$$A x = LU x = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ \frac{9}{25} & 1 & 0 & 0 \\ \frac{7}{10} & -\frac{45}{44} & 1 & 0 \\ \frac{3}{10} & -\frac{5}{44} & \frac{175}{57} & 1 \\ \end{array} \right)\left( \begin{array}{cccc} 50 & 60 & 23 & 60 \\ 0 & -\frac{88}{5} & -\frac{32}{25} & -\frac{38}{5} \\ 0 & 0 & \frac{57}{22} & -\frac{435}{22} \\ 0 & 0 & 0 & \frac{1935}{19} \\ \end{array} \right)x = y = \begin{pmatrix}40\\10\\67\\78\end{pmatrix}$$

The final result is

$$x = \begin{pmatrix}-\frac{1976}{645}\\-\frac{281}{860}\\ \frac{1327}{129}\\-\frac{256}{645} \end{pmatrix}$$

Related Solutions

[Math] Getting the inverse of a lower/upper triangular matrix

Ziyuang's answer handles the cases, where $N^2=0$, but it can be generalized as follows. A triangular $n\times n$ matrix $T$ with 1s on the diagonal can be written in the form $T=I+N$. Here $N$ is the strictly triangular part (with zeros on the diagonal), and it always satisfies the relation $N^{n}=0$. Therefore we can use the polynomial factorization $1-x^n=(1-x)(1+x+x^2+\cdots +x^{n-1})$ with $x=-N$ to get the matrix relation $$ (I+N)(I-N+N^2-N^3+\cdot+(-1)^{n-1}N^{n-1})=I + (-1)^{n-1}N^n=I $$ telling us that $(I+N)^{-1}=I+\sum_{k=1}^{n-1}(-1)^kN^k$.

Yet another way of looking at this is to notice that it also is an instance of a geometric series $1+q+q^2+q^3+\cdots =1/(1-q)$ with $q=-N$. The series converges for the unusual reason that powers of $q$ are all zero from some point on. The same formula can be used to good effect elsewhere in algebra, too. For example, in a residue class ring like $\mathbf{Z}/2^n\mathbf{Z}$ all the even numbers are nilpotent, so computing the modular inverse of an odd number can be done with this formula.

Linear Algebra – Solving Linear Matrix Equations of the Form $X = AXA^T + C$

Your solution for $X$ is correct.

Here is another solution (it is cheating, though, because I used Octave). On he plus side, it shows another approach to solving the equation.

$A$ is diagonalizable, and so has a basis of eigenvectors, $u_1,u_2$. Define the linear operator $L(X) = X-A X A^*$. You want to solve $L(X) = C$. Note that $u_i u_j^*$ forms a basis for the space of $2 \times 2 $ matrices, and $L(u_i u_j^*) = (1-\lambda_i \overline{\lambda_j}) u_i u_j^*$, hence $L$ is diagonalizable.

To express $C$ in terms of $u_i u_j^*$, let $C = \sum_{i,j} [\Gamma]_{ij} u_i u_j^*$, and solve for $\Gamma$. Let $U$ be the matrix whose columns are $u_1 ,u_2$. Note that $\sum_i [\Gamma]_{ij} u_i = U \Gamma e_j$, so we have $C = \sum_{j} U \Gamma e_j u_j^* = \sum_{j} U \Gamma e_j (U e_j)^* = \sum_{j} U \Gamma e_j e_j^T U^* = U \Gamma (\sum_{j} e_j e_j^T) U^* = U \Gamma U^*$. It follows that $\Gamma = U^{-1} C (U^{-1})^*$, and if we solve $L(\sum_{i,j} [\tilde{X}]_{ij} u_i u_j^*) = \sum_{i,j} [\Gamma]_{ij} u_i u_j^*$, we get $[\tilde{X}]_{ij} = \frac{1}{1-\lambda_i \overline{\lambda_j}} [\Gamma]_{ij}$.

Reversing the procedure to obtain $X$ in the standard basis, we have $X = U \tilde{X} U^*$.

Using Octave we obtain the same result as above.

a = [1.5 1 ; -0.7 0 ]
c = [1 0.5 ; 0.5 0.25 ]

[u,d] = eig(a)

l = [d(1,1) ; d(2,2) ]

phi = ones(2,2)-l*l'

gamma = inv(u)*c*inv(u')

x_tilde = (1./phi).*gamma
x = u*x_tilde*u'

# check result...
a*x*a'+c-x

Here is the output:

octave:44> a*x*a'+c-x
ans =

   0.00000 - 0.00000i   0.00000 - 0.00000i
   0.00000 - 0.00000i   0.00000 + 0.00000i

octave:45> x
x =

   18.8802 +  0.0000i  -11.3672 -  0.0000i
  -11.3672 -  0.0000i    9.5013 +  0.0000i

Best Answer

Related Solutions

[Math] Getting the inverse of a lower/upper triangular matrix

Linear Algebra – Solving Linear Matrix Equations of the Form $X = AXA^T + C$

Related Question