Solution of linear algebraic equations
PROBLEM
We wish to find a solution to the matrix equation A.*X.=B. where operator * denotes matrix multiplications, A. is a square matrix, B. is a known right-hand side vector, and X. is an unknown vector.
ALGORITHM
The
LUBKSB algorithm (for forward substitution and backsubstitution) solves the set of
N linear equations
A.*X.=B., here
A. is output from
LUDCMP routine. For given matrix
A. and vector
B. solves the following statements:
call LUDCMP N; call LUBKSB N
|
the matrix equation A.*X.=B. and the solution is saved in B. vector.
IMPLEMENTATION
Unit:
nonrecursive internal subroutine
Global variables: input N by N matrix A. - LU decomposition, input vector Indx. - permutation vector; A. and Indx. are determined by the LUDCMP routine
Parameters: a positive integer N
Result: the B. vector includes the solution of matrix equation given above
Note:
The A. matrix and Indx. vector are not modified. We can use them for successive calls with different B. vectors.
LUBKSB: procedure expose A. B. Indx.
parse arg N
L = 0
do I = 1 to N
P = Indx.I; Sum = B.P; B.P = B.I
if L <> 0
then do
do J = L to I - 1
Sum = Sum - A.I.J * B.J
end
end
else if Sum <> 0 then L = I
B.I = Sum
end
do I = N to 1 by -1
Sum = B.I
do J = I + 1 to N
Sum = Sum - A.I.J * B.J
end
B.I = Sum / A.I.I
end
return
|
EXAMPLE
N=3
A.1.1 = 2; A.1.2 = -1; A.1.3 = -1; B.1 = 4
A.2.1 = 3; A.2.2 = 4; A.2.3 = -2; B.2 = 11
A.3.1 = 3; A.3.2 = -2; A.3.3 = 4; B.3 = 11
call LUDCMP N; call LUBKSB N
S = ""
do J = 1 to N; S = S B.J; end
say S
|
displays
3 1 1
CONNECTIONS
Literature
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in C : the art of scientific computing
- 2nd ed. University Press, Cambridge, 1992
Faddejev A.K., Sominskij J.S. Sbornik zadac po vyssej algebre
Nauka, Moskva 1964
