LU computes the LU decomposition of a matrix, \(A\), such that \(P A = L U\),
where \(L\) is a lower triangle matrix, \(U\) is an upper triangle, and \(P\) is a
permutation matrix.
Usage
LU(A, b, tol = sqrt(.Machine$double.eps), verbose = FALSE, ...)Arguments
- A
coefficient matrix
- b
right-hand side vector. When supplied the returned object will also contain the solved \(d\) and
xelements- tol
tolerance for checking for 0 pivot
- verbose
logical; if
TRUE, print intermediate steps- ...
additional arguments passed to
showEqn
Value
A list of matrix components of the solution, P, L and U. If b
is supplied, the vectors \(d\) and x are also returned.
Details
The LU decomposition is used to solve the equation \(A x = b\) by calculating \(L(Ux - d) = 0\), where \(Ld = b\). If row exchanges are necessary for \(A\) then the permutation matrix \(P\) will be required to exchange the rows in \(A\); otherwise, \(P\) will be an identity matrix and the LU equation will be simplified to \(A = L U\).
Examples
A <- matrix(c(2, 1, -1,
-3, -1, 2,
-2, 1, 2), 3, 3, byrow=TRUE)
b <- c(8, -11, -3)
(ret <- LU(A)) # P is an identity; no row swapping
#> $P
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
#>
#> $L
#> [,1] [,2] [,3]
#> [1,] 1.0 0 0
#> [2,] -1.5 1 0
#> [3,] -1.0 4 1
#>
#> $U
#> [,1] [,2] [,3]
#> [1,] 2 1.0 -1.0
#> [2,] 0 0.5 0.5
#> [3,] 0 0.0 -1.0
#>
with(ret, L %*% U) # check that A = L * U
#> [,1] [,2] [,3]
#> [1,] 2 1 -1
#> [2,] -3 -1 2
#> [3,] -2 1 2
LU(A, b)
#> $P
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
#>
#> $L
#> [,1] [,2] [,3]
#> [1,] 1.0 0 0
#> [2,] -1.5 1 0
#> [3,] -1.0 4 1
#>
#> $U
#> [,1] [,2] [,3]
#> [1,] 2 1.0 -1.0
#> [2,] 0 0.5 0.5
#> [3,] 0 0.0 -1.0
#>
#> $d
#> [,1]
#> [1,] 8
#> [2,] 1
#> [3,] 1
#>
#> $x
#> [,1]
#> [1,] 2
#> [2,] 3
#> [3,] -1
#>
LU(A, b, verbose=TRUE)
#>
#> Initial equation:
#> 2*x1 + x2 - 1*x3 = 8
#> -3*x1 - 1*x2 + 2*x3 = -11
#> -2*x1 + x2 + 2*x3 = -3
#>
#> Lower triangle equation:
#> x1 = 8
#> -1.5*x1 + x2 = -11
#> -x1 + 4*x2 + x3 = -3
#>
#> Forward-solving operations:
#>
#> Equation: x1 = 8
#> Solution: x1 = 8
#>
#> Equation: -1.5*x1 + x2 = -11
#> Substitution: -1.5*8 + x2 = -11
#> Solution: x2 = (-11 - -12) = 1
#>
#> Equation: -x1 + 4*x2 + x3 = -3
#> Substitution: -8 + 4*1 + x3 = -3
#> Solution: x3 = (-3 - -4) = 1
#>
#> Intermediate solution: d = (8, 1, 1)
#>
#> Upper triangle equation:
#> 2*x1 + x2 - 1*x3 = 8
#> 0.5*x2 + 0.5*x3 = 1
#> - 1*x3 = 1
#>
#> Back-solving operations:
#>
#> Equation: - 1*x3 = 1
#> Solution: x3 = 1/-1 = -1
#>
#> Equation: 0.5*x2 + 0.5*x3 = 1
#> Substitution: 0.5*x2 + 0.5*-1 = 1
#> Solution: x2 = (1 - -0.5)/0.5 = 3
#>
#> Equation: 2*x1 + x2 - 1*x3 = 8
#> Substitution: 2*x1 + -1 - 1*3 = 8
#> Solution: x1 = (8 - 4)/2 = 2
#>
#> Final solution: x = (2, 3, -1)
LU(A, b, verbose=TRUE, fractions=TRUE)
#>
#> Initial equation:
#> 2*x1 + x2 - 1*x3 = 8
#> -3*x1 - 1*x2 + 2*x3 = -11
#> -2*x1 + x2 + 2*x3 = -3
#>
#> Lower triangle equation:
#> x1 = 8
#> -3/2*x1 + x2 = -11
#> -x1 + 4*x2 + x3 = -3
#>
#> Forward-solving operations:
#>
#> Equation: x1 = 8
#> Solution: x1 = 8
#>
#> Equation: -3/2*x1 + x2 = -11
#> Substitution: -3/2*8 + x2 = -11
#> Solution: x2 = (-11 - -12) = 1
#>
#> Equation: -x1 + 4*x2 + x3 = -3
#> Substitution: -8 + 4*1 + x3 = -3
#> Solution: x3 = (-3 - -4) = 1
#>
#> Intermediate solution: d = (8, 1, 1)
#>
#> Upper triangle equation:
#> 2*x1 + x2 - 1*x3 = 8
#> 1/2*x2 + 1/2*x3 = 1
#> - 1*x3 = 1
#>
#> Back-solving operations:
#>
#> Equation: - 1*x3 = 1
#> Solution: x3 = 1/-1 = -1
#>
#> Equation: 1/2*x2 + 1/2*x3 = 1
#> Substitution: 1/2*x2 + 1/2*-1 = 1
#> Solution: x2 = (1 - -1/2)/1/2 = 3
#>
#> Equation: 2*x1 + x2 - 1*x3 = 8
#> Substitution: 2*x1 + -1 - 1*3 = 8
#> Solution: x1 = (8 - 4)/2 = 2
#>
#> Final solution: x = (2, 3, -1)
# permutations required in this example
A <- matrix(c(1, 1, -1,
2, 2, 4,
1, -1, 1), 3, 3, byrow=TRUE)
b <- c(1, 2, 9)
(ret <- LU(A, b))
#> $P
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 0 1
#> [3,] 0 1 0
#>
#> $L
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 1 1 0
#> [3,] 2 0 1
#>
#> $U
#> [,1] [,2] [,3]
#> [1,] 1 1 -1
#> [2,] 0 -2 2
#> [3,] 0 0 6
#>
#> $d
#> [,1]
#> [1,] 1
#> [2,] 8
#> [3,] 0
#>
#> $x
#> [,1]
#> [1,] 5
#> [2,] -4
#> [3,] 0
#>
with(ret, P %*% A)
#> [,1] [,2] [,3]
#> [1,] 1 1 -1
#> [2,] 1 -1 1
#> [3,] 2 2 4
with(ret, L %*% U)
#> [,1] [,2] [,3]
#> [1,] 1 1 -1
#> [2,] 1 -1 1
#> [3,] 2 2 4