Lab 13: LU Cholesky Thomas Lab · Numerical Methods

Type [A | b] (or just A). Watch L fill with the multipliers while U is eliminated, then the forward and backward substitutions. Turn on pivoting to see PA = LU.

Cholesky works only for symmetric positive definite matrices. Try changing one diagonal entry until it fails.

LU once, many right-hand sides

Try this

Run it. Fill in the table from the opening slide.
Add a fourth reading: top edge 120°.
Print L and U and compare with the stepper.

Thomas algorithm

Try this

Run it: is the profile a straight line?
Set n = 9. Predict the answer first.
Add a heat source: d[i] += 10 for every i. How does the profile change?

Cholesky

Try this

Run it. Compare L with Problem 2.
Change A[2][2] to 1. What message do you get, and why?

Q1 Factor

\begin{pmatrix} 3 & 2 & 1 \\ 6 & 6 & 3 \\ 9 & 10 & 6\end{pmatrix}

by Doolittle and solve with

\mathbf{b} = (6, 15, 25)

.

L = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 2 & 1\end{pmatrix}

,

U = \begin{pmatrix} 3 & 2 & 1 \\ 0 & 2 & 1 \\ 0 & 0 & 1\end{pmatrix}

;

\mathbf{y} = (6, 3, 1)

;

\mathbf{x} = (1, 1, 1)

.

Q2 Why is LU better than Gaussian elimination when there are many right-hand sides?

Elimination costs

\tfrac23 n^3

and depends only on

A

; LU stores its result (L and U). Each further

\mathbf{b}

needs only forward and back substitution,

2n^2

operations, instead of repeating the

\tfrac23 n^3

.

Q3 When can Cholesky be used, and what does it cost?

When

A

is symmetric positive definite.

A = LL^T

costs about

\tfrac13 n^3

operations (half of LU) and needs no pivoting. If a number under a square root is

\le 0

,

A

is not positive definite.

Q4 Find the Cholesky factor of

\begin{pmatrix} 4 & -2 & 0 \\ -2 & 5 & -2 \\ 0 & -2 & 5\end{pmatrix}

.

l_{11} = 2

,

l_{21} = -1

,

l_{31} = 0

,

l_{22} = \sqrt{5 - 1} = 2

,

l_{32} = (-2 - 0)/2 = -1

,

l_{33} = \sqrt{5 - 1} = 2

.

Q5 Solve

2x_1 - x_2 = 1

,

-x_1 + 2x_2 - x_3 = 0

,

-x_2 + 2x_3 - x_4 = 0

,

-x_3 + 2x_4 = 1

by the Thomas algorithm.

c' = (-0.5, -0.6667, -0.75)

,

d' = (0.5, 0.3333, 0.25, 1)

; back:

x_4 = 1

,

x_3 = 0.25 + 0.75 = 1

,

x_2 = 1

,

x_1 = 1

.