delphi实现数学算法————LU分解法
procedure LUBKSB(var A:matrx2; N:integer; var INDX,B:array of real);
var
I,LL,J,II:integer;
SUM:real;
BEGIN
II:=0;
For I:=1 To N DO
begin
LL:=Trunc(INDX[I]);
Sum:=B[LL];
B[LL]:=B[I];
If II <> 0 Then
begin
For J:=II To I - 1 do
Sum:=Sum - A[I, J] * B[J];
end
Else If Sum <> 0 Then
II:=I;
B[I]:=Sum;
end;
For I:=N DownTo 1 do
begin
Sum:=B[I];
If I < N Then
begin
For J:=I + 1 To N do
Sum:=Sum - A[I, J] * B[J];
end;
B[I]:=Sum / A[I, I];
end;
end;
procedure LUDCMP(var A:matrx2;N:integer;var INDX:array of real;var D:real);
const
NMAX = 100; TINY = 1E-20;
var
VV:array[0..100] of real;
I,J,K,IMAX:Integer;
AAMAX,SUM,DUM:Real;
begin
D:=1;
For I:=1 To N do
begin
AAMAX:=0;
For J:=1 To N do
If (Abs(A[I, J]) > AAMAX) Then AAMAX:=Abs(A[I, J]);
If AAMAX = 0 Then
begin
ShowMessage('Singular matrix');
exit;
end;
VV[I]:= 1 / AAMAX;
end;
For J:=1 To N do
begin
If J > 1 Then
begin
For I:=1 To J - 1 do
begin
Sum:=A[I, J];
If I > 1 Then
begin
For K:=1 To I - 1 do
Sum:=Sum - A[I, K] * A[K, J];
A[I, J]:=Sum;
end;
end;
end;
AAMAX:=0;
For I:=J To N do
begin
Sum:=A[I, J];
If J > 1 Then
begin
For K:=1 To J - 1 do
Sum:=Sum - A[I, K] * A[K, J];
A[I, J]:=Sum;
end;
DUM:=VV[I] * Abs(Sum);
If DUM >= AAMAX Then
begin
IMAX:=I;
AAMAX:=DUM;
end;
end;
If J <> IMAX Then
begin
For K:=1 To N do
begin
DUM:=A[IMAX, K];
A[IMAX, K]:=A[J, K];
A[J, K]:=DUM;
end;
D:=-D;
VV[IMAX]:=VV[J];
end;
INDX[J]:=IMAX;
If J <> N Then
begin
If A[J, J] = 0 Then A[J, J]:=TINY;
DUM:=1 / A[J, J];
For I:=J + 1 To N do
A[I, J]:=A[I, J] * DUM;
end;
end;
If A[N, N]=0 Then A[N, N]:=TINY;
end;