192 lines
4.6 KiB
Plaintext
192 lines
4.6 KiB
Plaintext
Every square matrix <math>A</math> can be decomposed into a product of a lower triangular matrix <math>L</math> and a upper triangular matrix <math>U</math>,
|
|
as described in [[wp:LU decomposition|LU decomposition]].
|
|
|
|
:<math>A = LU</math>
|
|
|
|
It is a modified form of Gaussian elimination.
|
|
While the [[Cholesky decomposition]] only works for symmetric,
|
|
positive definite matrices, the more general LU decomposition
|
|
works for any square matrix.
|
|
|
|
There are several algorithms for calculating L and U.
|
|
To derive ''Crout's algorithm'' for a 3x3 example,
|
|
we have to solve the following system:
|
|
|
|
:<math>
|
|
A =
|
|
\begin{pmatrix}
|
|
a_{11} & a_{12} & a_{13}\\
|
|
a_{21} & a_{22} & a_{23}\\
|
|
a_{31} & a_{32} & a_{33}\\
|
|
\end{pmatrix}
|
|
=
|
|
\begin{pmatrix}
|
|
l_{11} & 0 & 0 \\
|
|
l_{21} & l_{22} & 0 \\
|
|
l_{31} & l_{32} & l_{33}\\
|
|
\end{pmatrix}
|
|
\begin{pmatrix}
|
|
u_{11} & u_{12} & u_{13} \\
|
|
0 & u_{22} & u_{23} \\
|
|
0 & 0 & u_{33}
|
|
\end{pmatrix}
|
|
= LU
|
|
</math>
|
|
|
|
We now would have to solve 9 equations with 12 unknowns. To make the system uniquely solvable, usually the diagonal elements of <math>L</math> are set to 1
|
|
|
|
:<math>l_{11}=1</math>
|
|
:<math>l_{22}=1</math>
|
|
:<math>l_{33}=1</math>
|
|
|
|
so we get a solvable system of 9 unknowns and 9 equations.
|
|
|
|
:<math>
|
|
A =
|
|
\begin{pmatrix}
|
|
a_{11} & a_{12} & a_{13}\\
|
|
a_{21} & a_{22} & a_{23}\\
|
|
a_{31} & a_{32} & a_{33}\\
|
|
\end{pmatrix}
|
|
=
|
|
\begin{pmatrix}
|
|
1 & 0 & 0 \\
|
|
l_{21} & 1 & 0 \\
|
|
l_{31} & l_{32} & 1\\
|
|
\end{pmatrix}
|
|
\begin{pmatrix}
|
|
u_{11} & u_{12} & u_{13} \\
|
|
0 & u_{22} & u_{23} \\
|
|
0 & 0 & u_{33}
|
|
\end{pmatrix}
|
|
=
|
|
\begin{pmatrix}
|
|
u_{11} & u_{12} & u_{13} \\
|
|
u_{11}l_{21} & u_{12}l_{21}+u_{22} & u_{13}l_{21}+u_{23} \\
|
|
u_{11}l_{31} & u_{12}l_{31}+u_{22}l_{32} & u_{13}l_{31} + u_{23}l_{32}+u_{33}
|
|
\end{pmatrix}
|
|
= LU
|
|
</math>
|
|
|
|
Solving for the other <math>l</math> and <math>u</math>, we get the following equations:
|
|
|
|
:<math>u_{11}=a_{11}</math>
|
|
:<math>u_{12}=a_{12}</math>
|
|
:<math>u_{13}=a_{13}</math>
|
|
|
|
:<math>u_{22}=a_{22} - u_{12}l_{21}</math>
|
|
:<math>u_{23}=a_{23} - u_{13}l_{21}</math>
|
|
|
|
:<math>u_{33}=a_{33} - (u_{13}l_{31} + u_{23}l_{32})</math>
|
|
|
|
and for <math>l</math>:
|
|
|
|
:<math>l_{21}=\frac{1}{u_{11}} a_{21}</math>
|
|
:<math>l_{31}=\frac{1}{u_{11}} a_{31}</math>
|
|
|
|
:<math>l_{32}=\frac{1}{u_{22}} (a_{32} - u_{12}l_{31})</math>
|
|
|
|
We see that there is a calculation pattern, which can be expressed as the following formulas, first for <math>U</math>
|
|
|
|
:<math>u_{ij} = a_{ij} - \sum_{k=1}^{i-1} u_{kj}l_{ik}</math>
|
|
|
|
and then for <math>L</math>
|
|
|
|
:<math>l_{ij} = \frac{1}{u_{jj}} (a_{ij} - \sum_{k=1}^{j-1} u_{kj}l_{ik})</math>
|
|
|
|
We see in the second formula that to get the <math>l_{ij}</math> below the diagonal, we have to divide by the diagonal element (pivot) <math>u_{jj}</math>, so we get problems when <math>u_{jj}</math> is either 0 or very small, which leads to numerical instability.
|
|
|
|
The solution to this problem is ''pivoting'' <math>A</math>, which means rearranging the rows of <math>A</math>, prior to the <math>LU</math> decomposition, in a way that the largest element of each column gets onto the diagonal of <math>A</math>. Rearranging the rows means to multiply <math>A</math> by a permutation matrix <math>P</math>:
|
|
|
|
:<math>PA \Rightarrow A'</math>
|
|
|
|
Example:
|
|
|
|
:<math>
|
|
\begin{pmatrix}
|
|
0 & 1 \\
|
|
1 & 0
|
|
\end{pmatrix}
|
|
\begin{pmatrix}
|
|
1 & 4 \\
|
|
2 & 3
|
|
\end{pmatrix}
|
|
\Rightarrow
|
|
\begin{pmatrix}
|
|
2 & 3 \\
|
|
1 & 4
|
|
\end{pmatrix}
|
|
</math>
|
|
|
|
The decomposition algorithm is then applied on the rearranged matrix so that
|
|
|
|
:<math>PA = LU</math>
|
|
|
|
|
|
;Task description:
|
|
The task is to implement a routine which will take a square nxn matrix <math>A</math> and return a lower triangular matrix <math>L</math>, a upper triangular matrix <math>U</math> and a permutation matrix <math>P</math>,
|
|
so that the above equation is fulfilled.
|
|
|
|
You should then test it on the following two examples and include your output.
|
|
|
|
|
|
;Example 1:
|
|
<pre>
|
|
A
|
|
|
|
1 3 5
|
|
2 4 7
|
|
1 1 0
|
|
|
|
L
|
|
|
|
1.00000 0.00000 0.00000
|
|
0.50000 1.00000 0.00000
|
|
0.50000 -1.00000 1.00000
|
|
|
|
U
|
|
|
|
2.00000 4.00000 7.00000
|
|
0.00000 1.00000 1.50000
|
|
0.00000 0.00000 -2.00000
|
|
|
|
P
|
|
|
|
0 1 0
|
|
1 0 0
|
|
0 0 1
|
|
</pre>
|
|
|
|
;Example 2:
|
|
<pre>
|
|
A
|
|
|
|
11 9 24 2
|
|
1 5 2 6
|
|
3 17 18 1
|
|
2 5 7 1
|
|
|
|
L
|
|
|
|
1.00000 0.00000 0.00000 0.00000
|
|
0.27273 1.00000 0.00000 0.00000
|
|
0.09091 0.28750 1.00000 0.00000
|
|
0.18182 0.23125 0.00360 1.00000
|
|
|
|
U
|
|
|
|
11.00000 9.00000 24.00000 2.00000
|
|
0.00000 14.54545 11.45455 0.45455
|
|
0.00000 0.00000 -3.47500 5.68750
|
|
0.00000 0.00000 0.00000 0.51079
|
|
|
|
P
|
|
|
|
1 0 0 0
|
|
0 0 1 0
|
|
0 1 0 0
|
|
0 0 0 1
|
|
</pre>
|
|
<br><br>
|
|
|