Computational Physics - Department of Physics

178 6 Linear Algebra

and since the determinat ofLis equal to 1 (by construction since the diagonals ofLequal 1)
we can use the inverse ofLto obtain

Ux=L−^1 w=y,

which yields the intermediate step
L−^1 w=y

and multiplying withLon both sides we reobtain Eq. (6.22). As soon as we haveywe can
obtainxthroughUx=y.
For our four-dimentional example this takes the form

y 1 =w 1

l 21 y 1 +y 2 =w 2

l 31 y 1 +l 32 y 2 +y 3 =w 3

l 41 y 1 +l 42 y 2 +l 43 y 3 +y 4 =w 4.

and

u 11 x 1 +u 12 x 2 +u 13 x 3 +u 14 x 4 =y 1

u 22 x 2 +u 23 x 3 +u 24 x 4 =y 2

u 33 x 3 +u 34 x 4 =y 3

u 44 x 4 =y 4

This example shows the basis for the algorithm needed to solve the set ofnlinear equations.
The algorithm goes as follows

• Set up the matrixAand the vectorwwith their correct dimensions. This determines
  the dimensionality of the unknown vectorx.


• Then LU decompose the matrixAthrough a call to the function

C++: ludcmp(double a, int n, int indx, double &d)

Fortran: CALL lu_decompose(a, n, indx, d)

This functions returns the LU decomposed matrixA, its determinant and the vector

indx which keeps track of the number of interchanges of rows.If the determinant is

zero, the solution is malconditioned.

• Thereafter you call the function

C++: lubksb(double a, int n, int indx, double w)

Fortran: CALL lu_linear_equation(a, n, indx, w)

which uses the LU decomposed matrixAand the vectorwand returnsxin the same

place asw. Upon exit the original content inwis destroyed. If you wish to keep this

information, you should make a backup of it in your calling function.

6.4.4 Inverse of a Matrix and the Determinant

The basic definition of the determinant ofAis