LU Decomposition function. Test for 5x5 is at the bottom.
LU_Decomposition <- function(A){
if(nrow(A)==2){
TwoSquareMatrixLU(A)
}
if(nrow(A)==3){
ThreeSquareMatrixLU(A)
}
if(nrow(A)==4){
FourSquareMatrixLU(A)
}
if(nrow(A)==5){
FiveSquareMatrixLU(A)
}
}
TwoSquareMatrixLU <- function(Matrixtwo){
p <- nrow(Matrixtwo)
multiplier1 <- (Matrixtwo[2,1]/Matrixtwo[1,1])
for(i in 1:3){
Matrixtwo[2,i] <- (-1*(Matrixtwo[1,i])*multiplier+Matrixtwo[2,i])
}
L = matrix(c(1,0,multiplier1,1), nrow = 2, ncol = 2, byrow=TRUE)
U = Matrixtwo
print("L is")
print(L)
print("U is")
print(U)
print(L %*% U)
return(Matrixthree)
}
ThreeSquareMatrixLU <- function(Matrixthree){
p <- nrow(Matrixthree)
cat("This matrix is a ", p, " by ", p, " square matrix" )
multiplier <- (Matrixthree[2,1]/Matrixthree[1,1])
for(i in 1:3){
Matrixthree[2,i] <- (-1*(Matrixthree[1,i])*multiplier+Matrixthree[2,i])
}
multiplier2 = (Matrixthree[3,1]/Matrixthree[1,1])
for(i in 1:3){
Matrixthree[3,i] <- (-1*(Matrixthree[1,i])*multiplier2+Matrixthree[3,i])
}
multiplier3 = (Matrixthree[3,2]/Matrixthree[2,2])
for(i in 2:3){
Matrixthree[3,i] <- (-1*(Matrixthree[2,i])*multiplier3+Matrixthree[3,i])
}
L = matrix(c(1,0,0,multiplier,1,0,multiplier2,multiplier3,1), nrow = 3, ncol = 3, byrow=TRUE)
U = Matrixthree
print("L is")
print(L)
print("U is")
print(U)
print(L %*% U)
return(Matrixthree)
}
FourSquareMatrixLU <- function(Matrixfour){
p <- nrow(Matrixfour)
print("ORIGINAL MATRIX")
print(Matrixfour)
## 2,1
multiplier1 <- ((Matrixfour[2,1]/Matrixfour[1,1]))
for(i in 1:4){
Matrixfour[2,i] <- (-1*(Matrixfour[1,i])*multiplier1+Matrixfour[2,i])
}
## 3,1
multiplier2 = (Matrixfour[3,1]/Matrixfour[1,1])
for(i in 1:4){
Matrixfour[3,i] <- (-1*(Matrixfour[1,i])*multiplier2+Matrixfour[3,i])
}
# 3,2
multiplier3 = (Matrixfour[3,2]/Matrixfour[2,2])
for(i in 2:4){
Matrixfour[3,i] <- (-1*(Matrixfour[2,i])*multiplier3+Matrixfour[3,i])
}
## 4,1
multiplier4 <- (Matrixfour[4,1]/Matrixfour[1,1])
for(i in 1:4){
Matrixfour[4,i] <- (-1*(Matrixfour[1,i])*multiplier4+Matrixfour[4,i])
}
## 4,2
multiplier5 = (Matrixfour[4,2]/Matrixfour[2,2])
for(i in 2:4){
Matrixfour[4,i] <- (-1*(Matrixfour[2,i])*multiplier5+Matrixfour[4,i])
}
# 4,3
multiplier6 = (Matrixfour[4,3]/Matrixfour[3,3])
for(i in 3:4){
Matrixfour[4,i] <- (-1*(Matrixfour[3,i])*multiplier6+Matrixfour[4,i])
}
L = matrix(c(1,0,0,0,multiplier1,1,0,0,multiplier2,multiplier3,1,0,multiplier4,multiplier5,multiplier6,1), nrow = 4, ncol = 4, byrow=TRUE)
U = Matrixfour
print("L is")
print(L)
print("U is")
print(U)
print("L times U is")
print(L %*% U)
return(Matrixfour)
}
FiveSquareMatrixLU <- function(Matrixfive){
p <- nrow(Matrixfive)
print("ORIGINAL MATRIX")
print(Matrixfive)
## 2,1
multiplier1 <- ((Matrixfive[2,1]/Matrixfive[1,1]))
for(i in 1:5){
Matrixfive[2,i] <- (-1*(Matrixfive[1,i])*multiplier1+Matrixfive[2,i])
}
## 3,1
multiplier2 = (Matrixfive[3,1]/Matrixfive[1,1])
for(i in 1:5){
Matrixfive[3,i] <- (-1*(Matrixfive[1,i])*multiplier2+Matrixfive[3,i])
}
# 3,2
multiplier3 = (Matrixfive[3,2]/Matrixfive[2,2])
for(i in 2:5){
Matrixfive[3,i] <- (-1*(Matrixfive[2,i])*multiplier3+Matrixfive[3,i])
}
## 4,1
multiplier4 <- (Matrixfive[4,1]/Matrixfive[1,1])
for(i in 1:5){
Matrixfive[4,i] <- (-1*(Matrixfive[1,i])*multiplier4+Matrixfive[4,i])
}
## 4,2
multiplier5 = (Matrixfive[4,2]/Matrixfive[2,2])
for(i in 2:5){
Matrixfive[4,i] <- (-1*(Matrixfive[2,i])*multiplier5+Matrixfive[4,i])
}
# 4,3
multiplier6 = (Matrixfive[4,3]/Matrixfive[3,3])
for(i in 3:5){
Matrixfive[4,i] <- (-1*(Matrixfive[3,i])*multiplier6+Matrixfive[4,i])
}
## 5,1
multiplier7 <- (Matrixfive[5,1]/Matrixfive[1,1])
for(i in 1:5){
Matrixfive[5,i] <- (-1*(Matrixfive[1,i])*multiplier7+Matrixfive[5,i])
}
## 5,2
multiplier8 = (Matrixfive[5,2]/Matrixfive[2,2])
for(i in 2:5){
Matrixfive[5,i] <- (-1*(Matrixfive[2,i])*multiplier8+Matrixfive[5,i])
}
# 5,3
multiplier9 = (Matrixfive[5,3]/Matrixfive[3,3])
for(i in 3:5){
Matrixfive[5,i] <- (-1*(Matrixfive[3,i])*multiplier9+Matrixfive[5,i])
}
# 5,4
multiplier10 = (Matrixfive[5,4]/Matrixfive[4,4])
for(i in 4:5){
Matrixfive[5,i] <- (-1*(Matrixfive[4,i])*multiplier10+Matrixfive[5,i])
}
L = matrix(c(1,0,0,0,0,multiplier1,1,0,0,0,multiplier2,multiplier3,1,0,0,multiplier4,multiplier5,multiplier6,1,0,multiplier7,multiplier8,multiplier9,multiplier10,1), nrow = 5, ncol = 5, byrow=TRUE)
U = Matrixfive
print("L is")
print(L)
print("U is")
print(U)
print("L times U is")
print(L %*% U)
}
5x5 Matrix test
testMatrixfive <- matrix(c(1,2,-3,4,5,4,1,0,3,0,0,2,1,2,-1,3,-3,0,5,2,2,1,2,3,-4),nrow = 5, ncol=5,byrow=TRUE)
LU_Decomposition(testMatrixfive)
## [1] "ORIGINAL MATRIX"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 2 -3 4 5
## [2,] 4 1 0 3 0
## [3,] 0 2 1 2 -1
## [4,] 3 -3 0 5 2
## [5,] 2 1 2 3 -4
## [1] "L is"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0.0000000 0.0000000 0.0000000 0
## [2,] 4 1.0000000 0.0000000 0.0000000 0
## [3,] 0 -0.2857143 1.0000000 0.0000000 0
## [4,] 3 1.2857143 -1.4516129 1.0000000 0
## [5,] 2 0.4285714 0.6451613 0.2321429 1
## [1] "U is"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 2 -3.000000e+00 4.000000 5.000000
## [2,] 0 -7 1.200000e+01 -13.000000 -20.000000
## [3,] 0 0 4.428571e+00 -1.714286 -6.714286
## [4,] 0 0 -8.881784e-16 7.225806 2.967742
## [5,] 0 0 0.000000e+00 0.000000 -1.785714
## [1] "L times U is"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 2 -3 4 5
## [2,] 4 1 0 3 0
## [3,] 0 2 1 2 -1
## [4,] 3 -3 0 5 2
## [5,] 2 1 2 3 -4