Simultaneous Linear Equations
Simultaneous Linear Equations
Simultaneous Linear Equations
Back Substitution:
Now the equations are solved starting from the last equation as it has only one unknown.
b ( n −1)
xn = n( n −1)
a nn
Then the second last equation, that is the ( n − 1) th equation, has two unknowns: xn and x n −1 ,
but xn is already known. This reduces the ( n − 1) th equation also to one unknown. Back
substitution hence can be represented for all equations by the formula
n
bi(i −1) − ∑ aij(i −1) x j
j =i +1
xi = for i = n − 1, n − 2,…,1
aii(i −1)
and
bn( n −1)
xn = ( n −1)
a nn
Example 1
The upward velocity of a rocket is given at three different times in Table 1.
Solution
Forward Elimination of Unknowns
Since there are three equations, there will be two steps of forward elimination of unknowns.
First step
Divide Row 1 by 25 and then multiply it by 64, that is, multiply Row 1 by 64/25 = 2.56 .
([25 5 1] [106.8])× 2.56 gives Row 1 as
[64 12.8 2.56] [273.408]
Subtract the result from Row 2
[64 8 1] [177.2]
− [64 12.8 2.56] [273.408]
0 − 4.8 − 1.56 − 96.208
to get the resulting equations as
⎡ 25 5 1 ⎤ ⎡ a1 ⎤ ⎡ 106.8 ⎤
⎢ 0 − 4.8 − 1.56⎥ ⎢a ⎥ = ⎢− 96.208⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣144 12 1 ⎥⎦ ⎢⎣ a3 ⎥⎦ ⎢⎣ 279.2 ⎥⎦
Divide Row 1 by 25 and then multiply it by 144, that is, multiply Row 1 by 144/25 = 5.76 .
([25 5 1] [106.8])× 5.76 gives Row 1 as
[144 28.8 5.76] [615.168]
Subtract the result from Row 3
[144 12 1] [279.2]
− [144 28.8 5.76 ] [615.168]
0 − 16.8 − 4.76 − 335.968
to get the resulting equations as
⎡25 5 1 ⎤ ⎡ a1 ⎤ ⎡ 106.8 ⎤
⎢ 0 − 4.8 − 1.56 ⎥ ⎢a ⎥ = ⎢ − 96.208 ⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 0 − 16.8 − 4.76⎥⎦ ⎢⎣ a3 ⎥⎦ ⎢⎣− 335.968⎥⎦
Second step
We now divide Row 2 by –4.8 and then multiply by –16.8, that is, multiply Row 2 by
− 16.8/ − 4.8 = 3.5 .
([0 − 4.8 −1.56] [− 96.208])× 3.5 gives Row 2 as
[0 − 16.8 − 5.46] [− 336.728]
Subtract the result from Row 3
[0 − 16.8 − 4.76] [− 335.968]
− [0 − 16.8 − 5.46] [− 336.728]
0 0 0.7 0.76
to get the resulting equations as
⎡25 5 1 ⎤ ⎡ a1 ⎤ ⎡ 106.8 ⎤
⎢ 0 − 4.8 − 1.56⎥ ⎢a ⎥ = ⎢− 96.208⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 0 0 0.7 ⎥⎦ ⎢⎣ a3 ⎥⎦ ⎢⎣ 0.76 ⎥⎦
Back substitution
From the third equation
0.7a3 = 0.76
0.76
a3 =
0.7
= 1.08571
Substituting the value of a3 in the second equation,
− 4.8a2 − 1.56a3 = −96.208
− 96.208 + 1.56 a3
a2 =
− 4.8
− 96.208 + 1.56 × 1.08571
=
− 4.8
= 19.6905
Substituting the value of a 2 and a3 in the first equation,
25a1 + 5a2 + a3 = 106.8
106.8 − 5a 2 − a3
a1 =
25
106.8 − 5 × 19.6905 − 1.08571
=
25
= 0.290472
Hence the solution vector is
⎡ a1 ⎤ ⎡0.290472⎤
⎢a ⎥ = ⎢ 19.6905 ⎥
⎢ 2⎥ ⎢ ⎥
⎢⎣ a3 ⎥⎦ ⎢⎣ 1.08571 ⎥⎦
The polynomial that passes through the three data points is then
v(t ) = a1t 2 + a 2 t + a3
= 0.290472t 2 + 19.6905t + 1.08571, 5 ≤ t ≤ 12
Since we want to find the velocity at t = 6, 7.5, 9 and 11 seconds, we could simply substitute
each value of t in v(t ) = 0.290472t 2 + 19.6905t + 1.08571 and find the corresponding
velocity. For example, at t = 6
2
v(6) = 0.290472(6) + 19.6905(6) + 1.08571
= 129.686 m/s
However we could also find all the needed values of velocity at t = 6, 7.5, 9, 11 seconds
using matrix multiplication.
⎡t 2 ⎤
⎢ ⎥
v(t ) = [0.290472 19.6905 1.08571] ⎢ t ⎥
⎢1⎥
⎣ ⎦
So if we want to find v(6), v(7.5), v(9), v(11), it is given by
⎡6 2 7.5 2 9 2 112 ⎤
[v(6) v(7.5) v(9) v(11)] = [0.290472 19.6905 1.08571] ⎢⎢ 6 7.5 9 11 ⎥
⎥
⎢1 1 1 1 ⎥⎦
⎣
⎡36 56.25 81 121⎤
= [0.290472 19.6905 1.08571] ⎢⎢ 6 7.5 9 11 ⎥⎥
⎢⎣ 1 1 1 1 ⎥⎦
= [129.686 165.104 201.828 252.828]
v(6) = 129.686 m/s
v(7.5) = 165.104 m/s
v(9) = 201.828 m/s
v(11) = 252.828 m/s
Example 2
Use Naïve Gauss elimination to solve
20 x1 + 15 x2 + 10 x3 = 45
− 3x1 − 2.249x2 + 7 x3 = 1.751
5x1 + x2 + 3x3 = 9
Use six significant digits with chopping in your calculations.
Solution
Working in the matrix form
⎡ 20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢− 3 − 2.249 7 ⎥ ⎢ x ⎥ = ⎢1.751⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 5 1 3 ⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣ 9 ⎥⎦
Forward Elimination of Unknowns
First step
Divide Row 1 by 20 and then multiply it by –3, that is, multiply Row 1 by − 3 / 20 = −0.15 .
([20 15 10] [45])× −0.15 gives Row 1 as
[− 3 − 2.25 −1.5] [− 6.75]
Subtract the result from Row 2
[− 3 − 2.249 7] [1.751]
− [− 3 − 2.25 − 1.5] [− 6.75]
0 0.001 8.5 8.501
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 0.001 8.5⎥ ⎢ x ⎥ = ⎢8.501⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 5 1 3 ⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣ 9 ⎥⎦
Divide Row 1 by 20 and then multiply it by 5, that is, multiply Row 1 by 5 / 20 = 0.25
([20 15 10] [45])× 0.25 gives Row 1 as
[5 3.75 2.5] [11.25]
Subtract the result from Row 3
[5 1 3] [9]
− [5 3.75 2.5] [11.25]
0 − 2.75 0.5 − 2.25
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 0.001 8.5⎥ ⎢ x ⎥ = ⎢ 8.501 ⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 0 − 2.75 0.5⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣− 2.25⎥⎦
Second step
Now for the second step of forward elimination, we will use Row 2 as the pivot equation and
eliminate Row 3: Column 2.
Divide Row 2 by 0.001 and then multiply it by –2.75, that is, multiply Row 2 by
− 2.75 / 0.001 = −2750 .
([0 0.001 8.5] [8.501])× −2750 gives Row 2 as
[0 − 2.75 − 23375] [− 23377.75]
Rewriting within 6 significant digits with chopping
[0 − 2.75 − 23375] [− 23377.7]
Subtract the result from Row 3
[0 − 2.75 0.5] [− 2.25]
− [0 − 2.75 − 23375] [− 23377.7]
0 0 23375.5 23375.45
Rewriting within 6 significant digits with chopping
[0 0 23375.5] [− 23375.4]
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 0.001 8.5 ⎥⎥ ⎢⎢ x 2 ⎥⎥ = ⎢⎢ 8.501 ⎥⎥
⎢
⎢⎣ 0 0 23375.5⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣23375.4⎥⎦
This is the end of the forward elimination steps.
Back substitution
We can now solve the above equations by back substitution. From the third equation,
23375.5x3 = 23375.4
23375.4
x3 =
23375.5
= 0.999995
Substituting the value of x3 in the second equation
0.001x2 + 8.5 x3 = 8.501
8.501 − 8.5 x3
x2 =
0.001
8.501 − 8.5 × 0.999995
=
0.001
8.501 − 8.49995
=
0.001
0.00105
=
0.001
= 1.05
Substituting the value of x3 and x 2 in the first equation,
20 x1 + 15 x2 + 10 x3 = 45
45 − 15 x 2 − 10 x3
x1 =
20
45 − 15 × 1.05 − 10 × 0.999995
=
20
45 − 15.75 − 9.99995
=
20
29.25 − 9.99995
=
20
19.2500
=
20
= 0.9625
Hence the solution is
⎡ x1 ⎤
[ X ] = ⎢⎢ x2 ⎥⎥
⎣⎢ x3 ⎦⎥
⎡ 0.9625 ⎤
= ⎢⎢ 1.05 ⎥⎥
⎢⎣0.999995⎥⎦
Round-off error: The Naïve Gauss elimination method is prone to round-off errors. This is
true when there are large numbers of equations as errors propagate. Also, if there is
subtraction of numbers from each other, it may create large errors. See the example below.
Example 3
Remember Example 2 where we used Naïve Gauss elimination to solve
20 x1 + 15 x2 + 10 x3 = 45
− 3x1 − 2.249x2 + 7 x3 = 1.751
5 x1 + x2 + 3x3 = 9
using six significant digits with chopping in your calculations? Repeat the problem, but now
use five significant digits with chopping in your calculations.
Solution
Writing in the matrix form
⎡ 20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢− 3 − 2.249 7 ⎥ ⎢ x ⎥ = ⎢1.751⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 5 1 3 ⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣ 9 ⎥⎦
Forward Elimination of Unknowns
First step
Divide Row 1 by 20 and then multiply it by –3, that is, multiply Row 1 by − 3 / 20 = −0.15 .
([20 15 10] [45])× −0.15 gives Row 1 as
[− 3 − 2.25 −1.5] [− 6.75]
Subtract the result from Row 2
[− 3 − 2.249 7] [1.751]
− [− 3 − 2.25 − 1.5] [− 6.75]
0 0.001 8.5 8.501
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 0.001 8.5⎥ ⎢ x ⎥ = ⎢8.501⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 5 1 3 ⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣ 9 ⎥⎦
Divide Row 1 by 20 and then multiply it by 5, that is, multiply Row 1 by 5 / 20 = 0.25 .
([20 15 10] [45])× 0.25 gives Row 1 as
[5 3.75 2.5] [11.25]
Subtract the result from Row 3
[5 1 3] [9]
− [5 3.75 2.5] [11.25]
0 − 2.75 0.5 − 2.25
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 0.001 8.5⎥ ⎢ x ⎥ = ⎢ 8.501 ⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 0 − 2.75 0.5⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣− 2.25⎥⎦
Second step
Now for the second step of forward elimination, we will use Row 2 as the pivot equation and
eliminate Row 3: Column 2.
Divide Row 2 by 0.001 and then multiply it by –2.75, that is, multiply Row 2 by
− 2.75 / 0.001 = −2750 .
([0 0.001 8.5] [8.501])× −2750 gives Row 2 as
[0 − 2.75 − 23375] [− 23377.75]
Rewriting within 5 significant digits with chopping
[0 − 2.75 − 23375] [− 23377]
Subtract the result from Row 3
[0 − 2.75 0.5] [− 2.25]
− [0 − 2.75 − 23375] [− 23377 ]
0 0 23375 23374
Rewriting within 6 significant digits with chopping
[0 0 23375] [− 23374]
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 0.001 8.5 ⎥ ⎢ x ⎥ = ⎢ 8.501 ⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 0 0 23375⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣23374⎥⎦
This is the end of the forward elimination steps.
Back substitution
We can now solve the above equations by back substitution. From the third equation,
23375 x3 = 23374
23374
x3 =
23375
= 0.99995
Substituting the value of x3 in the second equation
0.001x2 + 8.5 x3 = 8.501
8.501 − 8.5 x3
x2 =
0.001
8.501 − 8.5 × 0.99995
=
0.001
8.501 − 8.499575
=
0.001
8.501 − 8.4995
=
0.001
0.0015
=
0.001
= 1.5
Substituting the value of x3 and x 2 in the first equation,
20 x1 + 15 x2 + 10 x3 = 45
45 − 15 x 2 − 10 x3
x1 =
20
45 − 15 ×1.5 − 10 × 0.99995
=
20
45 − 22.5 − 9.9995
=
20
22.5 − 9.9995
=
20
12.5005
=
20
12.500
=
20
= 0.625
Hence the solution is
⎡ x1 ⎤
[X ] = ⎢⎢ x2 ⎥⎥
⎢⎣ x3 ⎥⎦
⎡ 0.625 ⎤
= ⎢⎢ 1.5 ⎥⎥
⎢⎣0.99995⎥⎦
Compare this with the exact solution of
⎡ x1 ⎤ ⎡1⎤
[X ] = ⎢⎢ x2 ⎥⎥ = ⎢⎢1⎥⎥
⎢⎣ x3 ⎥⎦ ⎢⎣1⎥⎦
What are some techniques for improving the Naïve Gauss elimination method?
As seen in Example 3, round off errors were large when five significant digits were used as
opposed to six significant digits. One method of decreasing the round-off error would be to
use more significant digits, that is, use double or quad precision for representing the
numbers. However, this would not avoid possible division by zero errors in the Naïve Gauss
elimination method. To avoid division by zero as well as reduce (not eliminate) round-off
error, Gaussian elimination with partial pivoting is the method of choice.
How does Gaussian elimination with partial pivoting differ from Naïve Gauss
elimination?
The two methods are the same, except in the beginning of each step of forward elimination, a
row switching is done based on the following criterion. If there are n equations, then there
are n − 1 forward elimination steps. At the beginning of the k th step of forward elimination,
one finds the maximum of
a kk , ak +1,k , …………, a nk
Then if the maximum of these values is a pk in the p th row, k ≤ p ≤ n , then switch rows p
and k .
The other steps of forward elimination are the same as the Naïve Gauss elimination method.
The back substitution steps stay exactly the same as the Naïve Gauss elimination method.
Example 4
In the previous two examples, we used Naïve Gauss elimination to solve
20 x1 + 15 x2 + 10 x3 = 45
− 3x1 − 2.249x2 + 7 x3 = 1.751
5 x1 + x2 + 3x3 = 9
using five and six significant digits with chopping in the calculations. Using five significant
digits with chopping, the solution found was
⎡ x1 ⎤
[X ] = ⎢⎢ x2 ⎥⎥
⎢⎣ x3 ⎥⎦
⎡ 0.625 ⎤
= ⎢⎢ 1.5 ⎥⎥
⎢⎣0.99995⎥⎦
This is different from the exact solution of
⎡ x1 ⎤
[X ] = ⎢⎢ x2 ⎥⎥
⎢⎣ x3 ⎥⎦
⎡1⎤
= ⎢⎢1⎥⎥
⎢⎣1⎥⎦
Find the solution using Gaussian elimination with partial pivoting using five significant digits
with chopping in your calculations.
Solution
⎡ 20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢− 3 − 2.249 7 ⎥ ⎢ x ⎥ = ⎢1.751⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 5 1 3 ⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣ 9 ⎥⎦
Forward Elimination of Unknowns
Now for the first step of forward elimination, the absolute value of the first column elements
below Row 1 is
20 , − 3 , 5
or
20, 3, 5
So the largest absolute value is in the Row 1. So as per Gaussian elimination with partial
pivoting, the switch is between Row 1 and Row 1 to give
⎡ 20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢− 3 − 2.249 7 ⎥ ⎢ x ⎥ = ⎢1.751⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 5 1 3 ⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣ 9 ⎥⎦
Divide Row 1 by 20 and then multiply it by –3, that is, multiply Row 1 by − 3 / 20 = −0.15 .
([20 15 10] [45])× −0.15 gives Row 1 as
[− 3 − 2.25 −1.5] [− 6.75]
Subtract the result from Row 2
[− 3 − 2.249 7] [1.751]
− [− 3 − 2.25 − 1.5] [− 6.75]
0 0.001 8.5 8.501
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 0.001 8.5⎥ ⎢ x ⎥ = ⎢8.501⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 5 1 3 ⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣ 9 ⎥⎦
Divide Row 1 by 20 and then multiply it by 5, that is, multiply Row 1 by 5 / 20 = 0.25 .
([20 15 10] [45])× 0.25 gives Row 1 as
[5 3.75 2.5] [11.25]
Subtract the result from Row 3
[5 1 3] [9]
− [5 3.75 2.5] [11.25]
0 − 2.75 0.5 − 2.25
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 0.001 8.5⎥ ⎢ x ⎥ = ⎢ 8.501 ⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 0 − 2.75 0.5⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣− 2.25⎥⎦
This is the end of the first step of forward elimination.
Now for the second step of forward elimination, the absolute value of the second column
elements below Row 1 is
0.001 , − 2.75
or
0.001, 2.75
So the largest absolute value is in Row 3. So Row 2 is switched with Row 3 to give
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 − 2.75 0.5⎥ ⎢ x ⎥ = ⎢− 2.25⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 0 0.001 8.5⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣ 8.501 ⎥⎦
Divide Row 2 by –2.75 and then multiply it by 0.001, that is, multiply Row 2 by
0.001 / − 2.75 = −0.00036363 .
([0 − 2.75 0.5] [− 2.25])× −0.00036363 gives Row 2 as
[0 0.00099998 − 0.00018182] [0.00081816]
Subtract the result from Row 3
[0 0.001 8.5] [8.501]
− [0 0.00099998 − 0.00018182] [0.00081816]
0 0 8.50018182 8.50018184
Rewriting within 5 significant digits with chopping
[0 0 8.5001] [8.5001]
to get the resulting equations as
⎡20 15 10 ⎤ ⎡ x1 ⎤ ⎡ 45 ⎤
⎢ 0 − 2.75 0.5 ⎥⎥ ⎢⎢ x 2 ⎥⎥ = ⎢⎢ − 2.25 ⎥⎥
⎢
⎢⎣ 0 0 8.5001⎥⎦ ⎢⎣ x3 ⎥⎦ ⎢⎣8.5001⎥⎦
Back substitution
8.5001x3 = 8.5001
8.5001
x3 =
8.5001
=1
Substituting the value of x3 in Row 2
− 2.75 x2 + 0.5 x3 = −2.25
− 2.25 − 0.5 x 2
x2 =
− 2.75
− 2.25 − 0.5 × 1
=
− 2.75
− 2.25 − 0.5
=
− 2.75
− 2.75
=
− 2.75
=1
Substituting the value of x3 and x 2 in Row 1
20 x1 + 15 x2 + 10 x3 = 45
45 − 15 x 2 − 10 x3
x1 =
20
45 − 15 × 1 − 10 × 1
=
20
45 − 15 − 10
=
20
30 − 10
=
20
20
=
20
=1
So the solution is
⎡ x1 ⎤
[X ] = ⎢⎢ x2 ⎥⎥
⎢⎣ x3 ⎥⎦
⎡1⎤
= ⎢1⎥
⎢⎥
⎢⎣1⎥⎦
This, in fact, is the exact solution. By coincidence only, in this case, the round-off error is
fully removed.
Can we use Naïve Gauss elimination methods to find the determinant of a square
matrix?
One of the more efficient ways to find the determinant of a square matrix is by taking
advantage of the following two theorems on a determinant of matrices coupled with Naïve
Gauss elimination.
Theorem 1:
Let [A] be a n × n matrix. Then, if [B] is a n × n matrix that results from adding or
subtracting a multiple of one row to another row, then det( A) = det( B) (The same is true for
column operations also).
Theorem 2:
Let [A] be a n × n matrix that is upper triangular, lower triangular or diagonal, then
det( A) = a11 × a22 × ... × aii × ... × ann
n
= ∏ aii
i =1
This implies that if we apply the forward elimination steps of the Naïve Gauss elimination
method, the determinant of the matrix stays the same according to Theorem 1. Then since at
the end of the forward elimination steps, the resulting matrix is upper triangular, the
determinant will be given by Theorem 2.
Example 5
Find the determinant of
⎡ 25 5 1⎤
[ A] = ⎢⎢ 64 8 1⎥⎥
⎢⎣144 12 1⎥⎦
Solution
Remember in Example 1, we conducted the steps of forward elimination of unknowns using
the Naïve Gauss elimination method on [A] to give
⎡25 5 1 ⎤
[B] = ⎢ 0 − 4.8 − 1.56⎥⎥
⎢
⎢⎣ 0 0 0.7 ⎥⎦
According to Theorem 2
det( A) = det( B)
= 25 × ( −4.8) × 0.7
= −84.00
What if I cannot find the determinant of the matrix using the Naïve Gauss elimination
method, for example, if I get division by zero problems during the Naïve Gauss
elimination method?
Well, you can apply Gaussian elimination with partial pivoting. However, the determinant of
the resulting upper triangular matrix may differ by a sign. The following theorem applies in
addition to the previous two to find the determinant of a square matrix.
Theorem 3:
Let [A] be a n × n matrix. Then, if [B] is a matrix that results from switching one row with
another row, then det( B) = − det( A) .
Example 6
Find the determinant of
⎡ 10 − 7 0⎤
[A] = ⎢− 3 2.099 6⎥⎥
⎢
⎢⎣ 5 − 1 5⎥⎦
Solution
The end of the forward elimination steps of Gaussian elimination with partial pivoting, we
would obtain
⎡10 − 7 0 ⎤
⎢
[B] = ⎢ 0 2.5 5 ⎥⎥
⎢⎣ 0 0 6.002⎥⎦
det(B) = 10 × 2.5 × 6.002
= 150.05
Since rows were switched once during the forward elimination steps of Gaussian elimination
with partial pivoting,
det(A) = − det(B)
= −150.05
Example 7
Prove
1
det( A) =
( )
det A−1
Solution
[ A][ A]−1 = [ I ]
( )
det A A −1 = det (I )
det ( A)det(A ) = 1
−1
1
det ( A) =
( )
det A−1
If [A] is a n × n matrix and det( A) ≠ 0 , what other statements are equivalent to it?
1. [A] is invertible.
2. [ A]−1 exists.
3. [ A][ X ] = [C ] has a unique solution.
!
4. [ A][ X ] = [0] solution is [ X ] = [0] .
5. [ A][ A]−1 = [ I ] = [ A]−1 [ A] .
Gauss-Seidel Method
Why do we need another method to solve a set of simultaneous linear equations?
In certain cases, such as when a system of equations is large, iterative methods of solving
equations are more advantageous. Elimination methods, such as Gaussian elimination, are
prone to large round-off errors for a large set of equations. Iterative methods, such as the
Gauss-Seidel method, give the user control of the round-off error. Also, if the physics of the
problem are well known, initial guesses needed in iterative methods can be made more
judiciously leading to faster convergence.
What is the algorithm for the Gauss-Seidel method? Given a general set of n equations and
n unknowns, we have
a11 x1 + a12 x2 + a13 x3 + ... + a1n xn = c1
a21 x1 + a22 x2 + a23 x3 + ... + a2 n xn = c2
. .
. .
. .
an1 x1 + an 2 x2 + an3 x3 + ... + ann xn = cn
If the diagonal elements are non-zero, each equation is rewritten for the corresponding
unknown, that is, the first equation is rewritten with x1 on the left hand side, the second
equation is rewritten with x 2 on the left hand side and so on as follows
Example 1
The upward velocity of a rocket is given at three different times in the following table
Iteration #2
The estimate of the solution vector at the end of Iteration #1 is
⎡ a1 ⎤ ⎡ 3.6720 ⎤
⎢a ⎥ = ⎢− 7.8510⎥
⎢ 2⎥ ⎢ ⎥
⎢⎣ a3 ⎥⎦ ⎢⎣ − 155.36 ⎥⎦
Now we get
106.8 − 5(− 7.8510) − (−155.36)
a1 =
25
= 12.056
177.2 − 64(12.056) − (−155.36)
a2 =
8
= −54.882
279.2 − 144(12.056) − 12(− 54.882)
a3 =
1
= − 798.34
The absolute relative approximate error for each xi then is
12.056 − 3.6720
∈a 1 = ×100
12.056
= 69.543%
− 54.882 − (− 7.8510)
∈a 2 = ×100
− 54.882
= 85.695%
− 798.34 − (− 155.36 )
∈a 3 = ×100
− 798.34
= 80.540%
At the end of the second iteration the estimate of the solution vector is
⎡ a1 ⎤ ⎡ 12.056 ⎤
⎢a ⎥ = ⎢− 54.882⎥
⎢ 2⎥ ⎢ ⎥
⎣⎢ a3 ⎦⎥ ⎣⎢− 798.54 ⎦⎥
and the maximum absolute relative approximate error is 85.695%.
Conducting more iterations gives the following values for the solution vector and the
corresponding absolute relative approximate errors.
Iteration a1 ∈a 1 % a2 ∈a 2 % a3 ∈a 3 %
1 3.6720 72.767 –7.8510 125.47 –155.36 103.22
2 12.056 69.543 –54.882 85.695 –798.34 80.540
3 47.182 74.447 –255.51 78.521 –3448.9 76.852
4 193.33 75.595 –1093.4 76.632 –14440 76.116
5 800.53 75.850 –4577.2 76.112 –60072 75.963
6 3322.6 75.906 –19049 75.972 –249580 75.931
As seen in the above table, the solution estimates are not converging to the true solution of
a1 = 0.29048
a2 = 19.690
a3 = 1.0857
The above system of equations does not seem to converge. Why?
Well, a pitfall of most iterative methods is that they may or may not converge. However, the
solution to a certain classes of systems of simultaneous equations does always converge
using the Gauss-Seidel method. This class of system of equations is where the coefficient
matrix [A] in [ A][ X ] = [C ] is diagonally dominant, that is
n
a ii ≥ ∑ a ij for all i
j =1
j ≠i
n
aii > ∑ aij for at least one i
j =1
j ≠i
If a system of equations has a coefficient matrix that is not diagonally dominant, it may or
may not converge. Fortunately, many physical systems that result in simultaneous linear
equations have a diagonally dominant coefficient matrix, which then assures convergence for
iterative methods such as the Gauss-Seidel method of solving simultaneous linear equations.
Example 2
Find the solution to the following system of equations using the Gauss-Seidel method.
12 x1 + 3x2 − 5x3 = 1
x1 + 5x2 + 3x3 = 28
3x1 + 7 x2 + 13x3 = 76
Use
⎡ x1 ⎤ ⎡1⎤
⎢ x ⎥ = ⎢0 ⎥
⎢ 2⎥ ⎢ ⎥
⎢⎣ x3 ⎥⎦ ⎢⎣1⎥⎦
as the initial guess and conduct two iterations.
Solution
The coefficient matrix
⎡12 3 − 5⎤
[A] = ⎢⎢ 1 5 3 ⎥⎥
⎢⎣ 3 7 13 ⎥⎦
is diagonally dominant as
a11 = 12 = 12 ≥ a12 + a13 = 3 + − 5 = 8
a22 = 5 = 5 ≥ a21 + a23 = 1 + 3 = 4
a33 = 13 = 13 ≥ a31 + a32 = 3 + 7 = 10
and the inequality is strictly greater than for at least one row. Hence, the solution should
converge using the Gauss-Seidel method.
Rewriting the equations, we get
1 − 3 x 2 + 5 x3
x1 =
12
28 − x1 − 3x3
x2 =
5
76 − 3 x1 − 7 x 2
x3 =
13
Assuming an initial guess of
⎡ x1 ⎤ ⎡1⎤
⎢ x ⎥ = ⎢0 ⎥
⎢ 2⎥ ⎢ ⎥
⎢⎣ x3 ⎥⎦ ⎢⎣1⎥⎦
Iteration #1
1 − 3(0) + 5(1)
x1 =
12
= 0.50000
28 − (0.50000 ) − 3(1)
x2 =
5
= 4.9000
76 − 3(0.50000 ) − 7(4.9000)
x3 =
13
= 3.0923
The absolute relative approximate error at the end of the first iteration is
0.50000 − 1
∈a 1 = ×100
0.50000
= 100.00%
4.9000 − 0
∈a 2 = ×100
4.9000
= 100.00%
3.0923 − 1
∈a 3 = × 100
3.0923
= 67.662%
The maximum absolute relative approximate error is 100.00%
Iteration #2
1 − 3(4.9000) + 5(3.0923)
x1 =
12
= 0.14679
28 − (0.14679 ) − 3(3.0923)
x2 =
5
= 3.7153
76 − 3(0.14679 ) − 7(3.7153)
x3 =
13
= 3.8118
At the end of second iteration, the absolute relative approximate error is
0.14679 − 0.50000
∈a 1 = ×100
0.14679
= 240.61%
3.7153 − 4.9000
∈a 2 = ×100
3.7153
= 31.889%
3.8118 − 3.0923
∈a 3
= × 100
3.8118
= 18.874%
The maximum absolute relative approximate error is 240.61%. This is greater than the value
of 100.00% we obtained in the first iteration. Is the solution diverging? No, as you conduct
more iterations, the solution converges as follows.
Iteration x1 ∈a 1 % x2 ∈a 2 % x3 ∈a 3 %
1 0.50000 100.00 4.9000 100.00 3.0923 67.662
2 0.14679 240.61 3.7153 31.889 3.8118 18.874
3 0.74275 80.236 3.1644 17.408 3.9708 4.0064
4 0.94675 21.546 3.0281 4.4996 3.9971 0.65772
5 0.99177 4.5391 3.0034 0.82499 4.0001 0.074383
6 0.99919 0.74307 3.0001 0.10856 4.0001 0.00101
Example 3
Given the system of equations
3x1 + 7 x2 + 13 x3 = 76
x1 + 5x2 + 3x3 = 28
12 x1 + 3x2 - 5x3 = 1
find the solution using the Gauss-Seidel method. Use
⎡ x1 ⎤ ⎡1⎤
⎢ x ⎥ = ⎢0 ⎥
⎢ 2⎥ ⎢ ⎥
⎢⎣ x3 ⎥⎦ ⎢⎣1⎥⎦
as the initial guess.
Solution
Rewriting the equations, we get
76 − 7 x2 − 13 x3
x1 =
3
28 − x1 − 3x3
x2 =
5
1 − 12 x1 − 3x2
x3 =
−5
Assuming an initial guess of
⎡ x1 ⎤ ⎡1⎤
⎢ x ⎥ = ⎢0 ⎥
⎢ 2⎥ ⎢ ⎥
⎢⎣ x3 ⎥⎦ ⎢⎣1⎥⎦
the next six iterative values are given in the table below.
Iteration x1 ∈a 1 % x2 ∈a 2 % x3 ∈a 3 %
1 21.000 95.238 0.80000 100.00 50.680 98.027
2 –196.15 110.71 14.421 94.453 –462.30 110.96
3 1995.0 109.83 –116.02 112.43 4718.1 109.80
4 –20149 109.90 1204.6 109.63 –47636 109.90
5 2.0364 × 105 109.89 –12140 109.92 4.8144 × 105 109.89
6 –2.0579 × 106 109.89 1.2272 × 105 109.89 –4.8653 × 106 109.89
You can see that this solution is not converging and the coefficient matrix is not diagonally
dominant. The coefficient matrix
⎡ 3 7 13 ⎤
[A] = ⎢⎢ 1 5 3 ⎥⎥
⎢⎣12 3 − 5⎥⎦
is not diagonally dominant as
a11 = 3 = 3 ≤ a12 + a13 = 7 + 13 = 20
Hence, the Gauss-Seidel method may or may not converge.
However, it is the same set of equations as the previous example and that converged. The
only difference is that we exchanged first and the third equation with each other and that
made the coefficient matrix not diagonally dominant.
Therefore, it is possible that a system of equations can be made diagonally dominant if one
exchanges the equations with each other. However, it is not possible for all cases. For
example, the following set of equations
x1 + x2 + x3 = 3
2 x1 + 3x2 + 4 x3 = 9
x1 + 7 x2 + x3 = 9
cannot be rewritten to make the coefficient matrix diagonally dominant.
LU Decomposition
I hear about LU decomposition used as a method to solve a set of simultaneous linear
equations. What is it?
We already studied two numerical methods of finding the solution to simultaneous linear
equations – Naïve Gauss elimination and Gaussian elimination with partial pivoting. Then,
why do we need to learn another method? To appreciate why LU decomposition could be a
better choice than the Gauss elimination techniques in some cases, let us discuss first what
LU decomposition is about.
For a nonsingular matrix [A] on which one can successfully conduct the Naïve Gauss
elimination forward elimination steps, one can always write it as
[A] = [L][U ]
where
[L]= Lower triangular matrix
[U ] = Upper triangular matrix
Then if one is solving a set of equations
[ A][ X ] = [C ] ,
then
[L][U ][X ] = [C] as ([ A] = [L][U ])
Multiplying both sides by [L ]−1 ,
[L]−1 [L][U ][X ] = [L]−1 [C ]
[I ][U ][X ]= [L]−1 [C ] as ([L]−1 [L] = [ I ])
[U ][X ] = [L]−1 [C ] as ([I ][U ] = [U ])
Let
[L]−1 [C ] = [Z ]
then
[L][Z ] = [C] (1)
and
[U ][X ] = [Z ] (2)
So we can solve Equation (1) first for [Z ] by using forward substitution and then use
Equation (2) to calculate the solution vector [X ] by back substitution.
This is all exciting but LU decomposition looks more complicated than Gaussian
elimination. Do we use LU decomposition because it is computationally more efficient
than Gaussian elimination to solve a set of n equations given by [A][X]=[C]?
For a square matrix [A] of n × n size, the computational time1 CT |DE to decompose the [A]
matrix to [ L][U ] form is given by
⎛ 8n 3 20n ⎞
CT |DE = T ⎜⎜ + 4n 2 − ⎟,
⎝ 3 3 ⎟⎠
where
T = clock cycle time2.
The computational time CT |FS to solve by forward substitution [L][Z ] = [C ] is given by
CT |FS = T (4n 2 − 4n )
The computational time CT |BS to solve by back substitution [U ][X ] = [Z ] is given by
CT |BS = T (4n 2 + 12n )
So, the total computational time to solve a set of equations by LU decomposition is
CT |LU = CT |DE + CT |FS + CT |BS
⎛ 8n 3 20n ⎞
= T ⎜⎜ + 4n 2 − ( ) (
⎟⎟ + T 4n 2 − 4n + T 4n 2 + 12n )
⎝ 3 3 ⎠
3
⎛ 8n 4n ⎞
= T ⎜⎜ + 12n 2 + ⎟⎟
⎝ 3 3⎠
Now let us look at the computational time taken by Gaussian elimination. The computational
time CT |FE for the forward elimination part,
⎛ 8n3 32n ⎞
CT |FE = T ⎜⎜ + 8n 2 − ⎟,
⎝ 3 3 ⎟⎠
and the computational time CT |BS for the back substitution part is
CT |BS = T (4n 2 + 12n )
So, the total computational time CT |GE to solve a set of equations by Gaussian Elimination
is
CT |GE = CT |FE + CT |BS
1
The time is calculated by first separately calculating the number of additions, subtractions,
http://www.isi.edu/~draper/papers/mwscas07_kwon.pdf
2
As an example, a 1.2 GHz CPU has a clock cycle of 1 /(1.2 ×109 ) = 0.833333 ns
⎛ 8n 3 32n ⎞
= T ⎜⎜ + 8n 2 − (
⎟⎟ + T 4n 2 + 12n )
⎝ 3 3 ⎠
3
⎛ 8n 4n ⎞
= T ⎜⎜ + 12n 2 + ⎟⎟
⎝ 3 3⎠
The computational time for Gaussian elimination and LU decomposition is identical.
This has confused me further! Why learn LU decomposition method when it takes the
same computational time as Gaussian elimination, and that too when the two methods
are closely related. Please convince me that LU decomposition has its place in solving
linear equations!
We have the knowledge now to convince you that LU decomposition method has its
place in the solution of simultaneous linear equations. Let us look at an example where the
LU decomposition method is computationally more efficient than Gaussian elimination.
Remember in trying to find the inverse of the matrix [A] in Chapter 04.05, the problem
reduces to solving n sets of equations with the n columns of the identity matrix as the RHS
vector. For calculations of each column of the inverse of the [A] matrix, the coefficient
matrix [A] matrix in the set of equation [A][X ] = [C] does not change. So if we use the LU
decomposition method, the [A] = [L][U ] decomposition needs to be done only once, the
forward substitution (Equation 1) n times, and the back substitution (Equation 2) n times.
Therefore, the total computational time CT |inverse LU required to find the inverse of a
matrix using LU decomposition is
CT |inverse LU = 1× CT |DE + n × CT |FS + n × CT |BS
⎛ 8n 3 20n ⎞
= 1 × T ⎜⎜ + 4n 2 − ( ) ( )
⎟⎟ + n × T 4n 2 − 4n + n × T 4n 2 + 12n
⎝ 3 3 ⎠
⎛ 32n 3 20n ⎞
= T ⎜⎜ + 12n 2 − ⎟⎟
⎝ 3 3 ⎠
In comparison, if Gaussian elimination method were used to find the inverse of a matrix, the
forward elimination as well as the back substitution will have to be done n times. The total
computational time CT |inverse GE required to find the inverse of a matrix by using Gaussian
elimination then is
CT |inverse GE = n × CT |FE + n × CT |BS
⎛ 8n 3 32n ⎞
= n × T ⎜⎜ + 8n 2 − (
⎟⎟ + n × T 4n 2 + 12n )
⎝ 3 3 ⎠
⎛ 8n 4 4n 2 ⎞
= T ⎜⎜ + 12n 3 + ⎟
⎝ 3 3 ⎟⎠
Clearly for large n , CT |inverse GE >> CT |inverse LU as CT |inverse GE has the dominating terms of n 4
and CT |inverse LU has the dominating terms of n 3 . For large values of n , Gaussian elimination
method would take more computational time (approximately n / 4 times – prove it) than the
LU decomposition method. Typical values of the ratio of the computational time for
different values of n are given in Table 1.
Are you convinced now that LU decomposition has its place in solving systems of equations?
We are now ready to answer other curious questions such as
1) How do I find LU matrices for a nonsingular matrix [A] ?
2) How do I conduct forward and back substitution steps of Equations (1) and (2),
respectively?
How do I decompose a non-singular matrix [A] , that is, how do I find [A] = [L][U]?
If forward elimination steps of the Naïve Gauss elimination methods can be applied on a
nonsingular matrix, then [A] can be decomposed into LU as
⎡ a11 … a1n ⎤
a12
⎢a ! a2n ⎥
a 22
[ A] = ⎢ 21 ⎥
⎢ " "
! " ⎥
⎢ ⎥
⎣ a n1 ! a nn ⎦
an 2
⎡1 0 … 0⎤ ⎡u11 u12 … u1n ⎤
⎢ℓ 1 ! 0⎥ ⎢ 0 u22 ! u2 n ⎥
⎢ 21 ⎥⎢ ⎥
=
⎢ " " ! "⎥ ⎢ " " ! " ⎥
⎢ ⎥⎢ ⎥
⎣ ℓ n1 ℓ n 2 ! 1 ⎦ ⎣ 0 0 ! unn ⎦
The elements of the [U ] matrix are exactly the same as the coefficient matrix one obtains at
the end of the forward elimination steps in Naïve Gauss elimination.
The lower triangular matrix [L ] has 1 in its diagonal entries. The non-zero elements on the
non-diagonal elements in [L ] are multipliers that made the corresponding entries zero in the
upper triangular matrix [U ] during forward elimination.
Let us look at this using the same example as used in Naïve Gaussian elimination.
Example 1
Find the LU decomposition of the matrix
⎡ 25 5 1⎤
[A] = ⎢⎢ 64 8 1⎥⎥
⎢⎣144 12 1⎥⎦
Solution
[A] = [L][U ]
⎡1 0 0⎤ ⎡u11 u12 u13 ⎤
= ℓ 21 1 0⎥ ⎢ 0 u22 u23 ⎥
⎢
⎢ ⎥⎢ ⎥
⎢⎣ℓ 31 ℓ 32 1⎥⎦ ⎢⎣ 0 0 u33 ⎥⎦
The [U ] matrix is the same as found at the end of the forward elimination of Naïve Gauss
elimination method, that is
⎡25 5 1 ⎤
[U ] = ⎢ 0 − 4.8 − 1.56⎥⎥
⎢
⎢⎣ 0 0 0.7 ⎥⎦
To find ℓ 21 and ℓ 31 , find the multiplier that was used to make the a 21 and a 31 elements zero
in the first step of forward elimination of the Naïve Gauss elimination method. It was
64
ℓ 21 =
25
= 2.56
144
ℓ 31 =
25
= 5.76
To find ℓ 32 , what multiplier was used to make a 32 element zero? Remember a 32 element
was made zero in the second step of forward elimination. The [A] matrix at the beginning of
the second step of forward elimination was
⎡25 5 1 ⎤
⎢ 0 − 4.8 − 1.56 ⎥
⎢ ⎥
⎢⎣ 0 − 16.8 − 4.76 ⎥⎦
So
− 16.8
ℓ 32 =
− 4.8
= 3.5
Hence
⎡ 1 0 0⎤
[L] = ⎢⎢2.56 1 0⎥⎥
⎢⎣5.76 3.5 1⎥⎦
Confirm [L][U ] = [A].
⎡ 1 0 0⎤ ⎡25 5 1 ⎤
[L][U ] = ⎢2.56 1 0⎥ ⎢ 0 − 4.8 − 1.56⎥⎥
⎢ ⎥ ⎢
⎢⎣5.76 3.5 1⎥⎦ ⎢⎣ 0 0 0.7 ⎥⎦
⎡ 25 5 1⎤
= ⎢⎢ 64 8 1⎥⎥
⎢⎣144 12 1⎥⎦
Example 2
Use the LU decomposition method to solve the following simultaneous linear equations.
⎡ 25 5 1⎤ ⎡ a1 ⎤ ⎡106.8 ⎤
⎢ 64 8 1⎥ ⎢a ⎥ = ⎢177.2 ⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣144 12 1⎥⎦ ⎢⎣ a3 ⎥⎦ ⎢⎣279.2⎥⎦
Solution
Recall that
[A][X ] = [C]
and if
[A] = [L][U ]
then first solving
[L][Z ] = [C]
and then
[U ][X ] = [Z ]
gives the solution vector [X ].
Now in the previous example, we showed
[A] = [L][U ]
⎡ 1 0 0⎤ ⎡25 5 1 ⎤
= ⎢2.56 1 0⎥ ⎢ 0 − 4.8 − 1.56⎥
⎢ ⎥⎢ ⎥
⎣⎢5.76 3.5 1⎦⎥ ⎣⎢ 0 0 0.7 ⎦⎥
First solve
[L][Z ] = [C]
⎡ 1 0 0⎤ ⎡ z1 ⎤ ⎡106.8 ⎤
⎢2.56 1 0⎥ ⎢ z ⎥ = ⎢177.2 ⎥
⎢ ⎥⎢ 2 ⎥ ⎢ ⎥
⎢⎣5.76 3.5 1⎥⎦ ⎢⎣ z 3 ⎥⎦ ⎢⎣279.2⎥⎦
to give
z1 = 106.8
2.56z1 + z2 = 177.2
5.76 z1 + 3.5z2 + z3 = 279.2
Forward substitution starting from the first equation gives
z1 = 106.8
z2 = 177.2 − 2.56z1
= 177.2 − 2.56 × 106.8
= −96.208
z3 = 279.2 − 5.76 z1 − 3.5z2
= 279.2 − 5.76 ×106.8 − 3.5 × (− 96.208)
= 0.76
Hence
⎡ z1 ⎤
[Z ] = ⎢⎢ z 2 ⎥⎥
⎢⎣ z 3 ⎥⎦
⎡ 106.8 ⎤
= ⎢⎢− 96.208⎥⎥
⎢⎣ 0.76 ⎥⎦
This matrix is same as the right hand side obtained at the end of the forward elimination steps
of Naïve Gauss elimination method. Is this a coincidence?
Now solve
[U ][X ] = [Z ]
⎡25 5 1 ⎤ ⎡ a1 ⎤ ⎡ 106.8 ⎤
⎢ 0 − 4.8 − 1.56 ⎥ ⎢a ⎥ = ⎢− 96.208⎥
⎢ ⎥ ⎢ 2⎥ ⎢ ⎥
⎢⎣ 0 0 0.7 ⎥⎦ ⎢⎣ a3 ⎥⎦ ⎢⎣ 0.76 ⎥⎦
25a1 + 5a2 + a3 = 106.8
− 4.8a2 − 1.56a3 = −96.208
0.7a3 = 0.76
From the third equation
0.7a3 = 0.76
0.76
a3 =
0.7
= 1.0857
Substituting the value of a3 in the second equation,
− 4.8a2 − 1.56a3 = −96.208
− 96.208 + 1.56 a3
a2 =
− 4.8
− 96.208 + 1.56 × 1.0857
=
− 4.8
= 19.691
Substituting the value of a 2 and a3 in the first equation,
25a1 + 5a2 + a3 = 106.8
106.8 − 5a 2 − a3
a1 =
25
106.8 − 5 × 19.691 − 1.0857
=
25
= 0.29048
Hence the solution vector is
⎡ a1 ⎤ ⎡0.29048⎤
⎢a ⎥ = ⎢ 19.691 ⎥
⎢ 2⎥ ⎢ ⎥
⎢⎣ a3 ⎥⎦ ⎢⎣ 1.0857 ⎥⎦
How do I find the inverse of a square matrix using LU decomposition?
A matrix [B] is the inverse of [A] if
[A][B] = [I ] = [B][A].
How can we use LU decomposition to find the inverse of the matrix? Assume the first
column of [B] (the inverse of [A]) is
[b11 b12 ... ... bn1 ]T
Then from the above definition of an inverse and the definition of matrix multiplication
⎡b11 ⎤ ⎡1⎤
⎢b ⎥ ⎢0⎥
[A]⎢ 21 ⎥ = ⎢ ⎥
⎢ ! ⎥ ⎢!⎥
⎢ ⎥ ⎢ ⎥
⎣bn1 ⎦ ⎣0⎦
Similarly the second column of [B] is given by
⎡ b12 ⎤ ⎡0⎤
⎢b ⎥ ⎢1⎥
[A]⎢ 22 ⎥ = ⎢ ⎥
⎢ ! ⎥ ⎢!⎥
⎢ ⎥ ⎢ ⎥
⎣bn 2 ⎦ ⎣0⎦
Similarly, all columns of [B] can be found by solving n different sets of equations with the
column of the right hand side being the n columns of the identity matrix.
Example 3
Use LU decomposition to find the inverse of
⎡ 25 5 1⎤
[A] = ⎢⎢ 64 8 1⎥⎥
⎢⎣144 12 1⎥⎦
Solution
Knowing that
[A] = [L][U ]
⎡ 1 0 0⎤ ⎡25 5 1 ⎤
= ⎢2.56 1 0⎥ ⎢ 0 − 4.8 − 1.56⎥⎥
⎢ ⎥ ⎢
⎢⎣5.76 3.5 1⎥⎦ ⎢⎣ 0 0 0.7 ⎥⎦
We can solve for the first column of [ B ] = [A]−1 by solving for
⎡ 25 5 1⎤ ⎡b11 ⎤ ⎡1⎤
⎢ 64 8 1⎥ ⎢b ⎥ = ⎢0⎥
⎢ ⎥ ⎢ 21 ⎥ ⎢ ⎥
⎢⎣144 12 1⎥⎦ ⎢⎣b31 ⎥⎦ ⎢⎣0⎥⎦
First solve
[L][Z ] = [C],
that is
⎡ 1 0 0⎤ ⎡ z1 ⎤ ⎡1⎤
⎢2.56 1 0⎥ ⎢ z ⎥ = ⎢0⎥
⎢ ⎥⎢ 2 ⎥ ⎢ ⎥
⎢⎣5.76 3.5 1⎥⎦ ⎢⎣ z 3 ⎥⎦ ⎢⎣0⎥⎦
to give
z1 = 1
2.56z1 + z 2 = 0
5.76 z1 + 3.5z 2 + z 3 = 0
Forward substitution starting from the first equation gives
z1 = 1
z2 = 0 − 2.56z1
= 0 − 2.56(1)
= −2.56
z 3 = 0 − 5.76 z1 − 3.5z 2
= 0 − 5.76(1) − 3.5(− 2.56)
= 3.2
Hence
⎡ z1 ⎤
[Z ] = ⎢⎢ z 2 ⎥⎥
⎢⎣ z 3 ⎥⎦
⎡ 1 ⎤
= ⎢⎢− 2.56 ⎥⎥
⎢⎣ 3.2 ⎥⎦
Now solve
[U ][X ] = [Z ]
that is
⎡25 5 1 ⎤ ⎡b11 ⎤ ⎡ 1 ⎤
⎢ 0 − 4.8 − 1.56 ⎥ ⎢b ⎥ = ⎢− 2.56⎥
⎢ ⎥ ⎢ 21 ⎥ ⎢ ⎥
⎢⎣ 0 0 0.7 ⎥⎦ ⎢⎣b31 ⎥⎦ ⎢⎣ 3.2 ⎥⎦
25b11 + 5b21 + b31 = 1
− 4.8b21 − 1.56b31 = −2.56
0.7b31 = 3.2
Backward substitution starting from the third equation gives
3.2
b31 =
0.7
= 4.571
− 2.56 + 1.56b31
b21 =
− 4 .8
− 2.56 + 1.56(4.571)
=
− 4.8
= −0.9524
1 − 5b21 − b31
b11 =
25
1 − 5(−0.9524) − 4.571
=
25
= 0.04762
Hence the first column of the inverse of [A] is
⎡b11 ⎤ ⎡ 0.04762 ⎤
⎢b ⎥ = ⎢− 0.9524⎥
⎢ 21 ⎥ ⎢ ⎥
⎢⎣b31 ⎥⎦ ⎢⎣ 4.571 ⎥⎦
Similarly by solving
⎡ 25 5 1⎤ ⎡b12 ⎤ ⎡0⎤ ⎡b12 ⎤ ⎡− 0.08333⎤
⎢ 64 8 1⎥ ⎢b ⎥ = ⎢1⎥ gives ⎢b ⎥ = ⎢ 1.417 ⎥
⎢ ⎥ ⎢ 22 ⎥ ⎢ ⎥ ⎢ 22 ⎥ ⎢ ⎥
⎢⎣144 12 1⎥⎦ ⎢⎣b32 ⎥⎦ ⎢⎣0⎥⎦ ⎢⎣b32 ⎥⎦ ⎢⎣ − 5.000 ⎥⎦
and solving
⎡ 25 5 1⎤ ⎡b13 ⎤ ⎡0⎤ ⎡b13 ⎤ ⎡ 0.03571 ⎤
⎢ 64 8 1⎥ ⎢b ⎥ = ⎢0⎥ gives ⎢b ⎥ = ⎢− 0.4643⎥
⎢ ⎥ ⎢ 23 ⎥ ⎢ ⎥ ⎢ 23 ⎥ ⎢ ⎥
⎢⎣144 12 1⎥⎦ ⎢⎣b33 ⎥⎦ ⎢⎣1⎥⎦ ⎢⎣b33 ⎥⎦ ⎢⎣ 1.429 ⎥⎦
Hence
⎡ 0.04762 − 0.08333 0.03571 ⎤
[A] = ⎢⎢− 0.9524 1.417 − 0.4643⎥⎥
−1