Ad
Table of contents

Why LU Decomposition Method.

As you know Gauss elimination is designed to solve systems of linear algebraic equations, [A]{X} = {B}

Although it certainly represents a sound way to solve such systems, it becomes inefficient when solving equations with the same coefficients [A], but with different constants ( b’s).

Ad

The Gauss elimination involves two steps: forward elimination and back substitution. Of these, the forward-elimination step comprises the bulk of the computational effort.

This is particularly true for large systems of equations. That is the reason why LU Decomposition method c++ has been designed.

Rule | LU Decomposition Method.

[A] {X} = {B}.
LU decomposition methods separate the time-consuming elimination of the matrix [A] from the manipulations of the right-hand side {B}.

Thus, once [A] has been “decomposed,” multiple right-hand-side vectors can be evaluated in an efficient manner.

LU Decomposition Method C++ Program

  //techindetail.com
//lu decomposition C++

#include<iostream>
#include<iomanip>

//size of matrix
# define N 3
typedef float matrix[N][N];
matrix l,u,a;

using namespace std;

float b[N],x[N],v[N];

//upper triangular matrix

void urow(int i)
{
	float s;
	int j,k;
	for(j=i;j<N;j++)
	{
		s=0;
		for(k=0;k<i;k++)
		  s+=u[k][j]*l[i][k];
		  u[i][j]=a[i][j]-s;
	}
}

//lower triangular matrix

void lcol(int j)
{
	float s;
	int i,k;
	for(i=j+1;i<N;i++)
	{
		s=0;
		for(k=0;k<=i-1;k++)
		s+=u[k][j]*l[i][k];
		l[i][j]=(a[i][j]-s)/u[j][j];
	}
}
//function to display the matrix

void printmat(matrix x)
{
	int i,j;
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		cout<<setw(8)<<setprecision(4)<<x[i][j];
		cout<<endl;
	}
}
int main()
{
	int i,j,m;
	float s;
	cout<<"Enter the elements of augmented matrix row wise"<<endl;
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		cin>>a[i][j];
		cin>>b[i];
	} 
	for(i=0;i<N;i++)
	l[i][i]=1.0;
	for(m=0;m<N;m++)
	{
		urow(m);
		if(m<(N-1))
		lcol(m);
	}
	cout<<setw(14)<<"U \n"<<endl;
	printmat(u);
	cout<<endl;
	cout<<setw(14)<<"L \n"<<endl;
	printmat(l);

	for(i=0;i<N;i++)
	{
		s=0;
		for(j=0;j<=i-1;j++)
		s+=l[i][j]*v[j];
		v[j]=b[i]-s;
	}
	for(i=N-1;i>=0;i--)
	{
		s=0;
		for(j=i+1;j<N;j++)
		s+=u[i][j]*x[j];
		x[i]=(v[i]-s)/u[i][i];
	}
	cout<<"The solution is: "<<endl;
	for(i=0;i<N;i++)
	cout<<"x["<<setw(3)<<i+1<<"] ="<<setw(6)<<setprecision(4)<<x[i]<<endl;
	return 0;
}Code language: C++ (cpp)

Mathematical Overview of LU Decomposition.

Just as was the case with Gauss elimination, LU decomposition requires pivoting to avoid division by zero. Let’s explain it using three simultaneous equations, then the result can be extended to n-dimensional system.
Let the equation be:
[A]{X} = {B}
The equation can be rearranged as:
[A]{X} − {B} = 0      (1)

The above equation can be expressed in upper triangular matrix as:
LU decomposition upper triangular matrix
The above matrix can be expressed in equation form as below:
[U]{X} = {D}    (2-1)
After rearranging the equation becomes
[U]{X} − {D} = 0      (2-2)

Now let’s assume that there is a lower triangular matrix with 1’s on diagonal.
i.e. LU decomposition lower triangular matrix

[L]{X} = {D}    (3)

=>[U]{X} – {D}=0    (4)
Now when we multiply the lower and upper triangular matrix equations we get the first equation back:
[L]{[U]{X} − {D}} = [A]{X} − {B}.

If this equation holds, it follows from the rules for matrix multiplication that,
[L][U] = [A]            5.
and
[L]{D} = {B}            6.

Ad

Two step strategy:

A two-step strategy for obtaining solutions can be based on Eqs.( 2-2),
(5), and (6):

  1. LU decomposition step. [A] is factored or “decomposed” into lower [L] and upper [U] triangular matrices.
  2. Substitution step. [L] and [U] are used to determine a solution {X} for a right-hand side {B}. This step itself consists of two steps.
  3. First, Eq. (6) is used to generate an intermediate vector {D} by forward substitution. Then, the result is substituted into Eq. (2-2), which can be solved by back substitution for {X}.

Some useful Methods