Introduction
The LU decomposition is the factorization of a square matrix into two triangular matrices, one upper triangular matrix, and one lower triangular matrix, with the product of these two matrices equaling the original matrix.
LU decomposition has a variety of applications, including solving a system of equations, which is an integral part of many other applications, such as finding current in a circuit, finding the inverse of a matrix, and finding the determinant of a matrix.
Doolittle’s Algorithm provides a way to factor a matrix A into an LU decomposition.
Doolittle Algorithm
A square matrix can always be factored into a lower triangular matrix and an upper triangular matrix. In other words, [A] = [L][U].
Doolittle's method allows us to factor A into a LU decomposition without going through the complexity of Gaussian Elimination.
We assume that a LU decomposition exists for a general n×n matrix A and write the form of L and U directly. We then systematically use the equations from the multiplications required for A = LU to solve the entries in L and U.
A = LU can be expanded as follow:
The terms of the U matrix is given by:
The terms of the L matrix is given by:
Example
Given a matrix A, decompose it to Lower and Upper triangular matrix using the Doolittle algorithm.
Solution:
Using Doolittle’s algorithm, we can write A = LU, where L is the lower triangular matrix and U is the upper triangular matrix.
Now, multiply the L and U matrix.
Now, compare the multiplied LU matrix with the A matrix row by row.
After comparing each row element in both the above matrices, we will get the values of unknowns as follows:
So finally, the L and U matrix becomes:
Code Example
The C++ program to calculate the Upper and Lower Triangular Matrix is given below:
#include <bits/stdc++.h>
using namespace std;
const int MAX = 100;
void LUDecomposition(int mat[][MAX], int n)
{
int lower[n][n], upper[n][n];
memset(lower, 0, sizeof(lower));
memset(upper, 0, sizeof(upper));
// Decomposing matrix into Upper and Lower triangular matrix
for (int i = 0; i < n; i++)
{
// Upper Triangular
for (int k = i; k < n; k++)
{
// Summation of L(i, j) * U(j, k)
int sum = 0;
for (int j = 0; j < i; j++)
sum += (lower[i][j] * upper[j][k]);
// Evaluating U(i, k)
upper[i][k] = mat[i][k] - sum;
}
// Lower Triangular
for (int k = i; k < n; k++)
{
if (i == k)
lower[i][i] = 1; // Diagonal as 1
else
{
// Summation of L(k, j) * U(j, i)
int sum = 0;
for (int j = 0; j < i; j++)
sum += (lower[k][j] * upper[j][i]);
// Evaluating L(k, i)
lower[k][i] = (mat[k][i] - sum) / upper[i][i];
}
}
}
// setw is for displaying in correct manner
cout << setw(6)
<< " Lower Triangular"
<< setw(32)
<< " Upper Triangular" << endl;
// Printing the result
for (int i = 0; i < n; i++)
{
// Lower Triangular Matrix
for (int j = 0; j < n; j++)
cout << setw(6) << lower[i][j] << "\t";
cout << "\t";
// Upper Triangular Matrix
for (int j = 0; j < n; j++)
cout << setw(6) << upper[i][j] << "\t";
cout << endl;
}
}
int main()
{
int mat[][MAX] = {{1, 1, 1}, {1, 2, 2}, {1, 2, 3}};
LUDecomposition(mat, 3);
return 0;
}
Output:
The matrix we decomposed is: