Основы компьютерного моделирования физических и технических систем. Основы моделирования физических и технических систем, страница 6

v:=x[i];

x[i]:=(b[i]-s)/a[i,i];

if abs(v-x[i])>m then

m:=abs(v-x[i]);

end;

inc(cnt);

if cnt = 10 then

Base := Log2(m/e)

else if (cnt > 10)and (cnt mod 20 = 0) then

begin

if M >= LastM then

exit;

SetProgress((Base-Log2(m/e))/Base);

LastM := M;

end;

until m<e;

end;

procedure gauss(N:Integer; A:TMatrix; B:TVector; var X:TVector);

var

h:Extended;

i,j,k:integer;

Begin {Начало основной программы}

{Прямой ход - исключение переменных}

for i:=1 to n-1 do

begin

for j:=i+1 to n do

begin

assert(a[i,i]<>0);

a[j,i]:=-a[j,i]/a[i,i];

for k:=i+1 to n do

a[j,k]:=a[j,k]+a[j,i]*a[i,k];

b[j]:=b[j]+a[j,i]*b[i]

end;

SetProgress(i/N);

end;

x[n]:=b[n]/a[n,n];

{Обратный ход - нахождение корней}

for i:=n-1 downto 1 do

begin

h:=b[i];

for j:=i+1 to n do

h:=h-x[j]*a[i,j];

x[i]:=h/a[i,i];

end;

End;

procedure SwapReal(var A, B: Extended);

var

R: Extended;

begin

R := A;

A := B;

B := R;

end;

procedure MyJordan(N:Integer; A:TMatrix; B:TVector; var X:TVector);

var

h:Extended;

i,j,k, BestRow:integer;

Begin {Начало основной программы}

{Прямой ход - исключение переменных}

for i:=1 to n do

begin

SetProgress(i/N);

{Выбор главного элемента по столбцу }

h := 0;

BestRow := -1;

for K := i to N do

if abs(A[K, I])>h then

begin

BestRow := K;

h := A[K, I];

end;

assert(BestRow>0);

if BestRow < 0 then

continue;

{Перемена местами уравнений}

if BestRow <> i then

begin

SwapReal(B[BestRow], B[I]);

for K := 1 to N do

SwapReal(A[BestRow, K], A[I, K]);

end;

if h = 0 then

continue;

for j:=1 to n do

if j <> i then

begin

h := -a[j, i]/a[i, i];

for k:=i+1 to n do

a[j,k]:=a[j,k]+h*a[i,k];

b[j]:=b[j]+h*b[i];

a[j,i]:=0;

end;

end;

for I := 1 to N do

if a[I,I] = 0 then

x[I] := -1

else

x[I]:=b[I]/a[I,I];

End;

procedure MultByReal(K: Extended; var R: TMatrix);

var

I, J, N, M: integer;

begin

GetSize(R, N, M);

for I := 1 to N do

for J := 1 to M do

R[I, J] := K*R[I, J];

end;

{ TSlau }

procedure TSlau.AddFullEquation(Koeffs: array of Extended; Value: Extended);

var

I: Integer;

begin

for I := 1 to Size do

A[CurEquation, I] := Koeffs[I];

B[CurEquation] := Value;

Inc(CurEquation);

end;

constructor TSlau.Create(aSize: Integer);

begin

inherited Create;

Size := aSize;

SetSize(Size, Size, A);

SetLength(B, Size+1);

CurEquation := 1;

end;

procedure TSlau.Solve(var X: TVector; aMethod: TMethod);

begin

SetLength(X, Size+1);

case aMethod of

mZeidel: zeidel(Size, 1e-6, A, B, X);

mGauss: gauss(Size, A, B, X);

mJordan: MyJordan(Size, A, B, X);

mLancosh: MyLancosh(Size, A, B, X);

end;

end;

procedure MultMatrixByVector(A: TMatrix; B: TVector; var R: TVector);

var

N, M, I, J: Integer;

begin

N := length(R)-1;

M := length(B)-1;

for I := 1 to N do

begin

R[I] := 0;

for J := 1 to M do

R[I] := R[I] + A[I,J]*B[J];

end;

end;

function ScalarProduct(V1, V2: TVector): Extended;

var

N, I: Integer;

begin

N := length(V1)-1;

Result := 0;

for I := 1 to N do

Result := Result + V1[I]*V2[I];

end;

procedure MyLancosh(N:Integer; A:TMatrix; B:TVector; var X:TVector);

var

AxS, S, R: TVector;

ai, bi, rp: Extended;

I,J: Integer;

Scan: Extended;

begin

{

x0=zeros(size(b));

x=x0;

r=A*x0-b;

s=A*x0-b;

while abs(s'*(A*s))>1e-9

rp=r'*r;

ai=rp/(s'*(A*s));

r=r-ai*A*s;

bi=(r'*r)/rp;

x=x-ai*s;

s=r+bi*s;

end

}

SetLength(S, N+1);

SetLength(R, N+1);

SetLength(AxS, N+1);

for I := 1 to N do

begin

R[i] := -B[I];

S[i] := -B[I];

X[I] := 0;

end;

MultMatrixByVector(A, S, AxS);

//  repeat

for J := 1 to N do

begin

rp := ScalarProduct(R, R);

ai := rp/ScalarProduct(S, AxS);

for I := 1 to N do

begin

X[I] := X[I]-ai*S[I];

R[I] := R[I]-ai*AxS[I];

end;

bi := ScalarProduct(R, R)/rp;

for I := 1 to N do

S[I] := R[I]+bi*S[I];

MultMatrixByVector(A, S, AxS);

Scan := ScalarProduct(S, AxS);

if Scan < 1e-9 then

break;

if J mod 10 = 0  then

SetProgress(J/N);

end;

end;

procedure InvGauss(N:Integer; const A:TMatrix; var InvA:TMatrix);

var

h:Extended;

i,j,k:integer;

Begin {Начало основной программы}

SetSize(N, N, InvA);

for I := 1 to N do

InvA[I, I] := 1;

{Прямой ход - приведение к верхнему треугольному виду}