Systems of Linear Algebraic Equations - Numerical Methods in Engineering with MATLAB - page 48

=

If i j (a diagonalterm) , the solutionis

A jj −

−

j 1

L jk , L jj = j = 2 , 3 ,..., n (2.19)

=

k 1

Foranondiagonaltermwe get

A i j −

L ik L jk

j − 1

L i j = / L jj , j = 2 , 3 ,..., n − 1, i = j + 1 , j + 2 ,..., n (2.20)

k = 1

choleski

Note that in Eqs. (2.19) and (2.20) A i j appearsonlyinthe formula for L i j . Therefore,
once L i j has been computed, A i j is nolongerneeded. This makes it possible to write
the elements of L over the lower triangular portion of A as theyarecomputed. The
elements above the principal diagonalof A will remain untouched. At the conclusion
of decomposition L isextractedwith the MATLABcommand tril(A) . If anegative L jj
isencountered during decomposition, an errormessage is printed and the program
isterminated.

functionL=choleski(A)
%ComputesLinCholeski'sdecompositionA=LL'.
%USAGE:L=choleski(A)

n = size(A,1);
forj=1:n
temp = A(j,j) - dot(A(j,1:j-1),A(j,1:j-1));
if temp < 0.0
error('Matrix is not positive definite')

end
A(j,j) = sqrt(temp);
fori=j+1:n
A(i,j)=(A(i,j) - dot(A(i,1:j-1),A(j,1:j-1)))/A(j,j);

end

end
L = tril(A)

Wecouldalso write the algorithm for forward and back substitutionsthat are
necessary in the solution of Ax

Butsince Choleski's decomposition has no ad-
vantages overDoolittle's decompositioninthe solution of simultaneousequations,
we will skip that part.

= b .