Solving Using LU Decomposition
Solving Using LU DecompositionLet's do an example.
In[113]:=
a={{2, 1, 0, 0}, {1, 2, 1, 0}, {0, 1, 2, 1}, {0, 0, 1, 2}}
Out[113]=
{{2, 1, 0, 0}, {1, 2, 1, 0}, {0, 1, 2, 1}, {0, 0, 1, 2}}
In[114]:=
MatrixForm[a]
Out[114]=
2 1 0 0 1 2 1 0 0 1 2 1 0 0 1 2
Let's solve Ax=b using LU decomposition. Mathematica will use back substitution to solve the system.
In[115]:=
b={2, 1, 4, 8}
Out[115]=
{2, 1, 4, 8}
What is going on here?
The function slu finds the LU Factorization without pivoting.
In[116]:=
?slu
slu[a_, l_, u_] yields a lower triangular matrix l and an
upper triangular matrix u such that a=lu. This is LU
factorization without pivoting.
In[117]:=
Clear[l,u]
slu[a,l,u]
MatrixForm[l]
MatrixForm[u]
Out[117]=
1 0 0 0
1
-
2 1 0 0
2
-
0 3 1 0
3
-
0 0 4 1
Out[118]=
2 1 0 0
3
-
0 2 1 0
4
-
0 0 3 1
5
-
0 0 0 4
In[119]:=
MatrixForm[a]
MatrixForm[l.u]
Out[119]=
2 1 0 0 1 2 1 0 0 1 2 1 0 0 1 2
Out[120]=
2 1 0 0 1 2 1 0 0 1 2 1 0 0 1 2
If we have to solve Ax=b, using the LU factorization changes the problem to LUx=b. Let c=Ux. Then we must solve Lc=b.
We are solving two equations, first solve Lc=b for c and then Ux=c for x. Why solve two equations when the original was only one? L and U are both very nice matrices and the two equations are easy to solve.
Let' s continue with the example.
In[121]:=
MatrixForm[l]
MatrixForm[b]
Out[121]=
1 0 0 0
1
-
2 1 0 0
2
-
0 3 1 0
3
-
0 0 4 1
Out[122]=
2 1 4 8
Note that Lc=b is easy to solve. This is called forward substitution, for reasons that should become clear.
c1=2
1/2 c1 + c2 = 1
2/3 c2 + c3 = 4
3/4 c3 + c4 = 8
Begin with c1=2, substitute into the next equation to get c2=0;
substitute into the next equation to get c3=4;
the next equation yields c4=5. So, c={2, 0, 4, 5}
In[123]:=
c=LinearSolve[l,b]
Out[123]=
{2, 0, 4, 5}
In[124]:=
MatrixForm[u]
MatrixForm[c]
Out[124]=
2 1 0 0
3
-
0 2 1 0
4
-
0 0 3 1
5
-
0 0 0 4
Out[125]=
2 0 4 5
Now solve Ux=c. This is called back substitution, for reasons that should become clear.
2x1 + x2 = 2
3/2 x2 + x3 = 0
4/3 x3 + x4 = 4
5/4 x4 =5
Begin with x4=4, substitute into the next equation to get x3=0;
substitute into the next equation to get x2=0;
the next equation yields x1=1. So, x={1, 0, 0, 4}
In[126]:=
x=LinearSolve[u,c]
Out[126]=
{1, 0, 0, 4}
Of course, as far as Mathematica is concerned we could have used:
In[127]:=
LinearSolve[a,b]
Out[127]=
{1, 0, 0, 4}
However, LU factorization has important theoretical and practical uses, as you will see.