Computational Physics - Department of Physics

206 6 Linear Algebra

This can be written in matrix form as

Ax=w.

We specialize here to the following case

−x 1 +x 2 − 4 x 3 = 0

2 x 1 + 2 x 2 = 1

3 x 1 + 3 x 2 + 2 x 3 =^12.

Obtain the solution (by hand) of this system of equations by doing Gaussian elimination.

2. Write therafter a program which implements Gaussian elimination (with pivoting) and
   solve the above system of linear equations. How many floatingpoint operations are in-
   volved in the solution via Gaussian elimination without pivoting? Can you estimate the
   number of floating point operations with pivoting?

6.2.If the matrixAis real, symmetric and positive definite, then it has a uniquefactorization
(called Cholesky factorization)
A=LU=LLT

whereLTis the upper matrix, implying that

LTi j=Lji.

The algorithm for the Cholesky decomposition is a special case of the general LU-decomposition
algorithm. The algorithm of this decomposition is as follows

• Calculate the diagonal elementLiiby setting up a loop fori= 0 toi=n− 1 (C++ indexing
  of matrices and vectors)

Lii=

(

Aii−

i− 1

∑

k= 0

L^2 ik

) 1 / 2

. (6.51)

• within the loop overi, introduce a new loop which goes fromj=i+ 1 ton− 1 and calculate

Lji=

1

Lii

(

Ai j−

i− 1

∑

k= 0

Likljk

)

. (6.52)

For the Cholesky algorithm we have always thatLii> 0 and the problem with exceedingly large
matrix elements does not appear and hence there is no need forpivoting. Write a function
which performs the Cholesky decomposition. Test your program against the standard LU
decomposition by using the matrix

A=





6 3 2

3 2 1

2 1 1



 (6.53)

Finally, use the Cholesky method to solve

0. 05 x 1 + 0. 07 x 2 + 0. 06 x 3 + 0. 05 x 4 = 0. 23

0. 07 x 1 + 0. 10 x 2 + 0. 08 x 3 + 0. 07 x 4 = 0. 32

0. 06 x 1 + 0. 08 x 2 + 0. 10 x 3 + 0. 09 x 4 = 0. 33

0. 05 x 1 + 0. 07 x 2 + 0. 09 x 3 + 0. 10 x 4 = 0. 31

You can also use the LU codes for linear equations to check theresults.

6.3.In this exercise we are going to solve the one-dimensional Poisson equation in terms of
linear equations.