MA 214: Introduction To Numerical Analysis: Shripad M. Garge. IIT Bombay (Shripad@math - Iitb.ac - In)
MA 214: Introduction To Numerical Analysis: Shripad M. Garge. IIT Bombay (Shripad@math - Iitb.ac - In)
MA 214: Introduction To Numerical Analysis: Shripad M. Garge. IIT Bombay (Shripad@math - Iitb.ac - In)
Lecture 38
Shripad M. Garge.
IIT Bombay
2021-2022
1/71
[email protected] MA214 (2021-2022) L38
System of linear equations
2/71
[email protected] MA214 (2021-2022) L38
Gaußian elimination method
We can then solve for the variable xn using the last equation.
Once the value of xn is found, we use it in the second-last equation
3/71
[email protected] MA214 (2021-2022) L38
Gaußian elimination method
The procedure will fail if the element aii is zero in the step of
eliminating xi because then the operation
ˆ ˙
aji
Ej ´ Ei Ñ pEj q
aii
The system may still have a solution, but the technique for finding
the solution must be altered. Let us understand this through an
example: ¨ ˛¨ ˛ ¨ ˛
1 ´1 2 x1 ´6
˝2 ´2 3‚˝x2 ‚ “ ˝´14‚.
1 1 1 x3 0
We perform pE2 ´ 2E1 q Ñ pE2 q and pE3 ´ E1 q Ñ pE3 q.
4/71
[email protected] MA214 (2021-2022) L38
Gaußian elimination method
The diagonal entry a22 is now zero. So the procedure can not
continue in its present form.
We search for the first nonzero entry below akk in the k-th column.
If apk “ 0 for each p, then the matrix A is not invertible and hence
the system does not have a unique solution. The procedure then
stops.
6/71
[email protected] MA214 (2021-2022) L38
Operation count
8/71
[email protected] MA214 (2021-2022) L38
Operation count
Mult/Div:
n´1
ÿ n´1
ÿ n´1
ÿ
pn ´ iqpn ´ i ` 2q “ pn ´ iq2 ` 2 pn ´ iq
i“1 i“1 i“1
n´1
ÿ n´1
ÿ
“ i2 ` 2 i
i“1 i“1
pn ´ 1qnp2n ´ 1q pn ´ 1qn
“ `2
6 2
2n3 ` 3n2 ´ 5n
“ .
6
This is the count of multiplication or division operations.
9/71
[email protected] MA214 (2021-2022) L38
Operation count
Add/Sub:
n´1
ÿ n´1
ÿ n´1
ÿ
pn ´ iqpn ´ i ` 1q “ pn ´ iq2 ` pn ´ iq
i“1 i“1 i“1
n´1
ÿ n´1
ÿ
“ i2 ` i
i“1 i“1
pn ´ 1qnp2n ´ 1q pn ´ 1qn
“ `
6 2
n3 ´ n
“ .
3
This is the count of addition or subtraction operations.
Mult/Div:
n´1
ÿ n´1
ÿ
` ˘
1` pn ´ iq ` 1 “ 1` pn ´ iq ` pn ´ 1q
i“1 i“1
n´1
ÿ n´1
ÿ
“ pn ´ iq ` n “ i `n
i“1 i“1
npn ´ 1q n2 ` n
“ `n “ .
2 2
11/71
[email protected] MA214 (2021-2022) L38
Operation count
Mult/Div:
2n3 ` 3n2 ´ 5n n2 ` n n3 n
` “ ` n2 ´
6 2 3 3
Add/Sub:
n3 ´ n n2 ´ n n3 n2 5n
` “ ` ´ .
3 2 3 2 6
12/71
[email protected] MA214 (2021-2022) L38
Operation count
n Mult/Div Add/Sub
3 17 11
10 430 375
50 44, 150 42, 875
100 3, 43, 400 3, 38, 250
13/71
[email protected] MA214 (2021-2022) L38
MA 214: Introduction to numerical analysis
Lecture 39
Shripad M. Garge.
IIT Bombay
2021-2022
14/71
[email protected] MA214 (2021-2022) L39
Pivoting
This row interchange has the form pEk q Ø pEp q, where p is the
smallest integer greater than k with apk ‰ 0.
15/71
[email protected] MA214 (2021-2022) L39
Pivoting
In our next example, we will see that even for small systems,
round-off error can dominate the calculations.
The disparity in the magnitudes of m21 a13 and a23 has introduced
a round-off error, but the round-off error has not yet been
propagated.
17/71
[email protected] MA214 (2021-2022) L39
Pivoting
x2 « 1.001
59.17 ´ p59.14qp1.001q
x1 « “ ´10
0.003
59.14
contains the small error of 0.001 multiplied by « 20000.
0.003
This ruins the approximation to the actual value x1 “ 10.
20/71
[email protected] MA214 (2021-2022) L39
Partial pivoting
24/71
[email protected] MA214 (2021-2022) L39
Scaled partial pivoting
Consequently,
The scaling factors are computed only once, so the second step
requires pn ´ 1q divisions and pn ´ 2q comparisons.
comparisons and
n
ÿ npn ` 1q 1
k“ ´ 1 “ pn ´ 1qpn ` 2q
k“2
2 2
27/71
[email protected] MA214 (2021-2022) L39
Operation count
Shripad M. Garge.
IIT Bombay
2021-2022
29/71
[email protected] MA214 (2021-2022) L40
LU decomposition
We now see that the steps used to solve a system of the form
Ax “ b can be used to factor the matrix A.
It is a matrix which differs from the identity matrix only in the first
column, below the diagonal element, where the entry 0 is replaced
by the negatives of multipliers, ´mj1 .
We then have
32/71
[email protected] MA214 (2021-2022) L40
LU decomposition
In general, with
Apkq x “ b pkq
already formed, we multiply by the matrix M pkq .
This is the matrix which differs from the identity matrix only in the
k-th column, below the diagonal entry, where 0 is replaced by
´mjk to obtain the next system
The process ends with the formation of Apnq x “ b pnq , where Apnq is
the upper triangular matrix, the result of the Gaußian elimination
method.
We have
Apnq “ M pn´1q M pn´2q ¨ ¨ ¨ M p1q A.
33/71
[email protected] MA214 (2021-2022) L40
LU decomposition
Note that each of the M pkq is lower triangular and so is the inverse
of each M pkq .
The matrix M pkq effects the row operation pEj ´ mjk Ek q Ñ pEj q
and its inverse Lpkq should reverse this operation, hence should
perform the operation pEj ` mjk Ek q Ñ pEj q.
Thus, the matrix Lpkq “ pM pkq q´1 differs from the identity matrix
only in the k-th column, below the diagonal entry, where 0 is
replaced by the multiplier, mjk , for j “ k ` 1, . . . , n.
34/71
[email protected] MA214 (2021-2022) L40
LU decomposition
35/71
[email protected] MA214 (2021-2022) L40
An example
We obtain the LU factorisation. The first operations are:
pE2 ´ 2E1 q Ñ pE2 q, pE3 ´ 3E1 q Ñ pE3 q, pE4 ´ p´1qE1 q Ñ pE4 q.
This gives
¨ ˛¨ ˛
1 0 0 0 1 1 0 3
˚2 1 0 0‹ ˚0 ´1 ´1 ´5‹ ` ˘` ˘
A“˚ ‹˚ ‹ “ Lp1q M p1q A .
˝3 0 1 0‚˝0 ´4 ´1 ´7‚
´1 0 0 1 0 3 3 2
Observe that we had the diagonal entry in the matrix L, lii “ 1 for
all i in this method.
38/71
[email protected] MA214 (2021-2022) L40
Permutation matrices
39/71
[email protected] MA214 (2021-2022) L40
Permutation matrices
Then we get
A “ P ´1 LU “ P 1 LU
which is a decomposition of A as a product of a permutation
matrix P 1 , a lower triangular matrix L and an upper triangular
matrix M, in that order.
41/71
[email protected] MA214 (2021-2022) L40
MA 214: Introduction to numerical analysis
Lecture 41
Shripad M. Garge.
IIT Bombay
2021-2022
42/71
[email protected] MA214 (2021-2022) L41
PLU decomposition
43/71
[email protected] MA214 (2021-2022) L41
PLU decomposition
1 2 0 2
0 1 1 0
1 0 0 0
45/71
[email protected] MA214 (2021-2022) L41
PLU decomposition
We have ¨ ˛
1 1 ´1 2
˚1 2 0 2‹
PA “ ˚
˝´1 ´1 2
‹.
0‚
0 0 ´1 1
This is the matrix on which we can perform the Gaußian
elimination, the operations being the same as the ones above:
46/71
[email protected] MA214 (2021-2022) L41
PLU decomposition
¨ ˛¨ ˛
1 0 0 0 1 1 ´1 2
˚1 1 0 0‹ ˚0
‹ ˚ 1 1 0‹
PA “ ˚ ‹ “ LU.
˝´1 0 1 0‚˝0 0 1 2‚
0 0 ´1 1 0 0 0 3
The final decomposition for A as a product of a permutation
matrix, a lower triangular matrix and an upper triangular matrix
can be obtained as A “ P ´1 LU:
¨ ˛¨ ˛¨ ˛
0 0 0 1 1 0 0 0 1 1 ´1 2
˚1 0 0 1‹ ˚ 1 1 0 0‹ ˚0 1 1 0‹
˝0 0 1 0‚˝´1 0 1 0‚˝0 0 1 2‚.
˚ ‹˚ ‹˚ ‹
0 1 0 0 0 0 ´1 1 0 0 0 3
48/71
[email protected] MA214 (2021-2022) L41
Diagonally dominant matrices
49/71
[email protected] MA214 (2021-2022) L41
Positive definite matrices
The matrix ¨ ˛
2 ´1 0
A “ ˝´1 2 ´1‚
0 ´1 2
is a positive definite matrix.
50/71
[email protected] MA214 (2021-2022) L41
Positive definite matrices
Further,
x t Ax ą 0
A is invertible,
aii ą 0 for each i,
paij q2 ă aii ajj for each i ‰ j,
max |akj | ď max |aii |.
1ďk,jďn 1ďiďn
for some 1 ď k ď n.
55/71
[email protected] MA214 (2021-2022) L41
Cholesky factorisation
If
¨ ˛ ¨ ˛¨ ˛
4 ´1 1 l11 0 0 l11 l21 l31
A “ ˝´1 4.25 2.75‚ “ ˝l21 l22 0 ‚˝ 0 l22 l32 ‚
1 2.75 3.5 l31 l32 l33 0 0 l33
then we get
¨ ˛ ¨ 2 ˛
4 ´1 1 l11 l11 l21 l11 l31
˝´1 4.25 2.75‚ “ ˝l11 l21 2 ` l2
l21 l21 l31 ` l22 l32 ‚
22
1 2.75 3.5 2 ` l2 ` l2
l11 l31 l21 l31 ` l22 l32 l31 32 33
56/71
[email protected] MA214 (2021-2022) L41
Cholesky factorisation
a11 : 2
4 “ l11 ùñ l11 “ 2,
a21 : ´1 “ l11 l21 ùñ l21 “ ´0.5,
a31 : 1 “ l11 l31 ùñ l31 “ 0.5,
a22 2 ` l2
: 4.25 “ l21 ùñ l22 “2
22
a32 : 2.75 “ l21 l31 ` l22 l32 ùñ l32 “ 1.5
a33 2 ` l2 ` l2
: 3.5 “ l31 ùñ l33 “1
32 33
and hence
¨ ˛ ¨ ˛¨ ˛
4 ´1 1 2 0 0 2 ´0.5 0.5
A “ ˝´1 4.25 2.75‚ “ ˝´0.5 2 0‚˝0 2 1.5‚ “ LLt .
1 2.75 3.5 0.5 1.5 1 0 0 1
57/71
[email protected] MA214 (2021-2022) L41
MA 214: Introduction to numerical analysis
Lecture 42
Shripad M. Garge.
IIT Bombay
2021-2022
58/71
[email protected] MA214 (2021-2022) L42
Iterative methods
59/71
[email protected] MA214 (2021-2022) L42
Distance in Rn
The l2 -distance:
˜ ¸1{2
n
ÿ
2
}x ´ y}2 “ pxi ´ yi q .
i“1
The l8 -distance:
which has the exact solution x “ p1, 1, 1qt but the Gaußian
elimination method performed using five-digit rounding arithmetic
and partial pivoting produces the approximate solution
61/71
[email protected] MA214 (2021-2022) L42
Distance in Rn
We have
˘1{2
p0.2001q2 ` p0.00009q2 ` p0.07462q2
`
}x ´ x̃}2 “
“ 0.21356
and
}x ´ x̃}8 “ maxt0.2001, 0.00009, 0.07462u
“ 0.2001.
Here, the second and the third components are good
approximations but the first component is a poor approximation
and that error dominates the computations of both the above
distances.
We now define:
63/71
[email protected] MA214 (2021-2022) L42
Norms on matrices
We have
3
ÿ 3
ÿ
|a1j | “ 1 ` 2 ` 1 “ 4, |a2j | “ 0 ` 3 ` 1 “ 4
j“1 j“1
and
3
ÿ
|a3j | “ 5 ` 1 ` 1 “ 7.
j“1
Hence }A}8 “ 7.
64/71
[email protected] MA214 (2021-2022) L42
Eigenvalues, eigenvectors
65/71
[email protected] MA214 (2021-2022) L42
Eigenvalues, eigenvectors
66/71
[email protected] MA214 (2021-2022) L42
Eigenvalues, eigenvectors
λ “ 3: pA ´ 3I qpvq “ 0 gives
¨ ˛¨ ˛ ¨ ˛
´1 0 0 v1 0
˝ 1 ´2 2‚˝v2 ‚ “ ˝0‚.
1 ´1 1 v3 0
This gives
v1 “ 0 and v2 “ v3
hence an eigenvector for eigenvalue 3 is of the form cp0, 1, 1qt for
0 ‰ c P R.
Note that even though our A has real entries, some roots of the
characteristic polynomial of A may be (non-real) complex numbers.
If A is an n ˆ n matrix then
“ ‰1{2
}A}2 “ ρpAt Aq
ρpAq ď }A}2 and
ρpAq ď }A}8 .
Let ¨ ˛
1 1 0
A “ ˝ 1 2 1‚.
´1 1 2
We apply the above result, and for that we need to compute
ρpAt Aq.
Here
¨ ˛¨ ˛ ¨ ˛
1 1 ´1 1 1 0 3 2 ´1
At A “ ˝1 2 1 ‚˝ 1 2 1‚ “ ˝ 2 6 4 ‚
0 1 2 ´1 1 2 ´1 4 5
?
and the eigenvalues of At A are 0, 7 ˘ 7. Hence
b
a ?
t
}A}2 “ ρpA Aq “ 7 ` 7 « 3.10576.
69/71
[email protected] MA214 (2021-2022) L42
Convergent matrices
For instance, if
ˆ ˙ ˆ ˙
1{2 0 k 1{2k 0
A“ then A “ .
1{4 1{2 k{2k`1 1{2k
70/71
[email protected] MA214 (2021-2022) L42
Convergent matrices
1 A is a convergent matrix.
2 limnÑ8 }An }2 “ 0.
3 limnÑ8 }An }8 “ 0.
4 ρpAq ă 1.
5 limnÑ8 An v “ 0 for every vector v.
71/71
[email protected] MA214 (2021-2022) L42