## 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).

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**

```
//wikkihut.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;
}
```

## 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: 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)

[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.

**[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.

## Two step strategy:

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

(5), and (6):

- LU decomposition step. [A] is factored or “decomposed” into lower [L] and upper [U] triangular matrices.
- 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.
- 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

**Gauss Elimination with Partial Pivoting****Gauss Jordon Method**- Regula Falsi Method
- Bisection Method
- Gauss Elimination Method C++