array_to_LaTeX = function(arr){
rows = apply(arr, MARGIN=1, paste, collapse = " & ")
matrix_string <- paste(rows, collapse = " \\\\ ")
return(paste("\\begin{bmatrix}", matrix_string, "\\end{bmatrix}"))
}Esempio di fattorizzazione LU
Appunti di calcolo numerico
Il metodo di eliminazione di Gauss, applicato ad un sistema lineare della forma:
\left\{ \begin{matrix} x_2 + 3 x_3 + 4 x_4 &= 6 \\ 2x_1 + 2x_2 + 4 x_4 &= 10 \\ x_2 + 3 x_3 + x_4 &= 3 \\ x_1 + 3 x_3 - x_4 &= 0 \\ \end{matrix} \right.
riscritto nella forma \bm{A}\bm{x} = \bm{b}
\bm{A}= \begin{pmatrix} 0 & 1 & 3 & 4 \\ 2 & 2 & 0 & 4 \\ 0 & 1 & 3 & 1 \\ 1 & 0 & 3 & -1 \end{pmatrix}, \quad \bm{x} = \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix}, \quad \bm{b} = \begin{pmatrix} 6 \\ 10 \\ 3 \\ 0 \end{pmatrix},
Alloco il sistema lineare in R
# R program to create a matrix
A = matrix(
# Taking sequence of elements
c(0, 1, 3, 4,
2, 2, 0, 4,
0, 1, 3, 1,
1, 0, 3, -1),
nrow = 4, # No of rows
ncol = 4, # No of columns
byrow = TRUE # arrange elements in row-wise order
)
b = matrix( c(6, 10, 3, 0), nrow = 4, ncol = 1, byrow = TRUE )
P = matrix( c(1,2,3,4), nrow = 4, ncol = 1, byrow = TRUE )
cat("The 4x4 matrix:\n")The 4x4 matrix:print(A) [,1] [,2] [,3] [,4]
[1,] 0 1 3 4
[2,] 2 2 0 4
[3,] 0 1 3 1
[4,] 1 0 3 -1conosco gia la soluzione, controllo il residuo:
xsol = matrix( c(1, 2, 0, 1), nrow = 4, ncol = 1, byrow = TRUE )
rownames(xsol) = c("x_1", "x_2", "x_3", "x_4")
res = b - A %*% xsol;
print(res) [,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0Inizio l’algoritmo di Gauss
Step 1 scambio righe
S1 = matrix(
c(0, 1, 0, 0,
1, 0, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1),
nrow=4,ncol=4,byrow=TRUE
)
A1 = S1 %*% A;
P1 = S1 %*% P;
print(A1) [,1] [,2] [,3] [,4]
[1,] 2 2 0 4
[2,] 0 1 3 4
[3,] 0 1 3 1
[4,] 1 0 3 -1print(P1) [,1]
[1,] 2
[2,] 1
[3,] 3
[4,] 4\bm{S}_1\bm{A} = \begin{bmatrix} 2 & 2 & 0 & 4 \\ 0 & 1 & 3 & 4 \\ 0 & 1 & 3 & 1 \\ 1 & 0 & 3 & -1 \end{bmatrix}
Faccio lo stesso lavoro sulla tabella TABLE
TABLE = A
TABLE = S1 %*% TABLE;
print(TABLE) [,1] [,2] [,3] [,4]
[1,] 2 2 0 4
[2,] 0 1 3 4
[3,] 0 1 3 1
[4,] 1 0 3 -1Step 1b applica eliminazione
v1 = matrix( c(0, 0, 0, 1/2), nrow=4,ncol=1, byrow=TRUE )
e1T = matrix( c(1, 0, 0, 0), nrow=1,ncol=4, byrow=TRUE )
ID = matrix(
c(1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1),
nrow=4,ncol=4,byrow=TRUE
)
L1 = ID - v1 %*% e1T;
A1b = L1 %*% A1;
print(L1) [,1] [,2] [,3] [,4]
[1,] 1.0 0 0 0
[2,] 0.0 1 0 0
[3,] 0.0 0 1 0
[4,] -0.5 0 0 1print(A1b) [,1] [,2] [,3] [,4]
[1,] 2 2 0 4
[2,] 0 1 3 4
[3,] 0 1 3 1
[4,] 0 -1 3 -3\bm{A}_1=\bm{L}_1\bm{S}_1\bm{A} = \begin{bmatrix} 2 & 2 & 0 & 4 \\ 0 & 1 & 3 & 4 \\ 0 & 1 & 3 & 1 \\ 0 & -1 & 3 & -3 \end{bmatrix}
Applico eliminazione sulla tabella TABLE
TABLE[2:4,1] = TABLE[2:4,1] / TABLE[1,1]
TABLE[2,2:4] = TABLE[2,2:4]-TABLE[2,1]*TABLE[1,2:4]
TABLE[3,2:4] = TABLE[3,2:4]-TABLE[3,1]*TABLE[1,2:4]
TABLE[4,2:4] = TABLE[4,2:4]-TABLE[4,1]*TABLE[1,2:4]
print(TABLE) [,1] [,2] [,3] [,4]
[1,] 2.0 2 0 4
[2,] 0.0 1 3 4
[3,] 0.0 1 3 1
[4,] 0.5 -1 3 -3Step 2 nessuno scambio
L’elemento A_{2,2}\neq 0 quindi non sevre newssuno scambio di righe.
S2 = ID
P2 = S2 %*% P1;
A2 = S2 %*% A1b;
print(P2) [,1]
[1,] 2
[2,] 1
[3,] 3
[4,] 4Step 2b applica eliminazione
v2 = matrix( c(0, 0, 1, -1), nrow=4,ncol=1, byrow=TRUE )
e2T = matrix( c(0, 1, 0, 0), nrow=1,ncol=4, byrow=TRUE )
L2 = ID - v2 %*% e2T;
A2b = L2 %*% A2;
print(L2) [,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 -1 1 0
[4,] 0 1 0 1print(A2b) [,1] [,2] [,3] [,4]
[1,] 2 2 0 4
[2,] 0 1 3 4
[3,] 0 0 0 -3
[4,] 0 0 6 1\bm{A}_2=\bm{L}_2\bm{S}_2\bm{L}_1\bm{S}_1\bm{A} = \begin{bmatrix} 2 & 2 & 0 & 4 \\ 0 & 1 & 3 & 4 \\ 0 & 0 & 0 & -3 \\ 0 & 0 & 6 & 1 \end{bmatrix}
Applico eliminazione sulla tabella TABLE
# Attenzione con byrow=TRUE indici riga e colonna sono scambiati
TABLE[3:4,2] = TABLE[3:4,2] / TABLE[2,2]
TABLE[3,3:4] = TABLE[3,3:4]-TABLE[3,2]*TABLE[2,3:4]
TABLE[4,3:4] = TABLE[4,3:4]-TABLE[4,2]*TABLE[2,3:4]
print(TABLE) [,1] [,2] [,3] [,4]
[1,] 2.0 2 0 4
[2,] 0.0 1 3 4
[3,] 0.0 1 0 -3
[4,] 0.5 -1 6 1Step 3 scambio righe
S3 = matrix(
c(1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 0, 1,
0, 0, 1, 0),
nrow=4,ncol=4,byrow=TRUE
)
A3 = S3 %*% A2b;
P3 = S3 %*% P2;
print(A3) [,1] [,2] [,3] [,4]
[1,] 2 2 0 4
[2,] 0 1 3 4
[3,] 0 0 6 1
[4,] 0 0 0 -3print(P3) [,1]
[1,] 2
[2,] 1
[3,] 4
[4,] 3\bm{A}_3=\bm{S}_3\bm{L}_2\bm{S}_2\bm{L}_1\bm{S}_1\bm{A} = \begin{bmatrix} 2 & 2 & 0 & 4 \\ 0 & 1 & 3 & 4 \\ 0 & 0 & 6 & 1 \\ 0 & 0 & 0 & -3 \end{bmatrix}
Applico scambio sulla tabella TABLE
TABLE = S3 %*% TABLE
print(TABLE) [,1] [,2] [,3] [,4]
[1,] 2.0 2 0 4
[2,] 0.0 1 3 4
[3,] 0.5 -1 6 1
[4,] 0.0 1 0 -3Step 3b eliminazione
Matrice hia in forma triangolare \bm{L}_3 = \bm{I} non serve fare niente.
L3 = ID;
A3b = L3 %*% A3;
print(L3) [,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1\bm{A}_3=\bm{L}_3\bm{S}_3\bm{L}_2\bm{S}_2\bm{L}_1\bm{S}_1\bm{A} = \begin{bmatrix} 2 & 2 & 0 & 4 \\ 0 & 1 & 3 & 4 \\ 0 & 0 & 6 & 1 \\ 0 & 0 & 0 & -3 \end{bmatrix}
Applico eliminazione sulla tabella TABLE
TABLE[4,3] = TABLE[4,3] / TABLE[3,3]
TABLE[4,4] = TABLE[4,4]-TABLE[4,3]*TABLE[3,4]
print(TABLE) [,1] [,2] [,3] [,4]
[1,] 2.0 2 0 4
[2,] 0.0 1 3 4
[3,] 0.5 -1 6 1
[4,] 0.0 1 0 -3Abbiamo quindi ottenuto
\begin{aligned}&\underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}}_{\bm{L}_3}\underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix}}_{\bm{S}_3}\underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 1 & 0 & 1 \end{bmatrix}}_{\bm{L}_2}\underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}}_{\bm{S}_2}\\&\underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ -0.5 & 0 & 0 & 1 \end{bmatrix}}_{\bm{L}_1}\underbrace{\begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}}_{\bm{S}_1}\underbrace{\begin{bmatrix} 0 & 1 & 3 & 4 \\ 2 & 2 & 0 & 4 \\ 0 & 1 & 3 & 1 \\ 1 & 0 & 3 & -1 \end{bmatrix}}_{\bm{A}}=\underbrace{\begin{bmatrix} 2 & 2 & 0 & 4 \\ 0 & 1 & 3 & 4 \\ 0 & 0 & 6 & 1 \\ 0 & 0 & 0 & -3 \end{bmatrix}}_{\bm{A}_3}\end{aligned}
LL = TABLE;
LL[upper.tri(TABLE)] = 0;
diag(LL) = 1;
UU = TABLE;
UU[lower.tri(TABLE)] = 0;
print(LL) [,1] [,2] [,3] [,4]
[1,] 1.0 0 0 0
[2,] 0.0 1 0 0
[3,] 0.5 -1 1 0
[4,] 0.0 1 0 1print(UU) [,1] [,2] [,3] [,4]
[1,] 2 2 0 4
[2,] 0 1 3 4
[3,] 0 0 6 1
[4,] 0 0 0 -3print(LL %*% UU) [,1] [,2] [,3] [,4]
[1,] 2 2 0 4
[2,] 0 1 3 4
[3,] 1 0 3 -1
[4,] 0 1 3 1