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

begin

Result := Points[I].X;

end;

function TTria.Y(I: Integer): Extended;

begin

Result := Points[I].Y;

end;

{ TRealPoint }

constructor TRealPoint.Create(aID: integer; ax, ay: Extended);

begin

inherited Create;

ID := aID;

X := ax;

Y := aY;

DX := 0;

DY := 0;

end;

end.

//Модуль: операции с матрицами и решение СЛАУ

unit uSLAU;

interface

uses Math, Classes;

type

//Перечисление методов решения слау

TMethod = (mZeidel, mGauss, mJordan, mLancosh);

//Тип - квадратная матрица произвольного размера

TMatrix = array of array of Extended;

//Тип - вектор произвольной длины

TVector = array of Extended;

//Класс - СЛАУ

TSlau = class

//Матрица А

A: TMatrix;

//Вектор B

B: TVector;

//Номер текущего добавляемого уравнения

CurEquation: Integer;

//Размер системы

Size: Integer;

//Конструктор - создает систему заданного размера

constructor Create(aSize: Integer);

//Добавляем уравнение с заданными коэффициентами и правой частью в систему

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

//Решает СЛАУ выбранным методом, сохраняя ответ в X

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

end;

//Процедура решения методом Зейделя

procedure zeidel(N:Integer; e:Extended; A:TMatrix; B:TVector; var X:TVector);

//Процедура решения методом Гаусса

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

//Процедура решения методом Гаусса-Жордана с выбором главного элемента по столбцу

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

//Процедура решения методом Ланцоша

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

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

//Перемножает матрицы A и B и сохраняет результат в R

procedure MultMatrix(A, B: TMatrix; var R: TMatrix);

//Транспонирует матрицу А и сохраняет результат в R

procedure Transpose(A: TMatrix; var R: TMatrix);

//Устанавливает размер матрицы R

procedure SetSize(N, M: Integer; var R: TMatrix);

//Возвращает размер матрицы R

procedure GetSize(R: TMatrix; var N, M: Integer);

//Умножает матрицу на число

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

//Меняет местами значения двух переменных

procedure SwapReal(var A, B: Extended);

//Умножает матрицу A на вектор B

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

//Возвращает скалярное произведение двух векторов

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

implementation

uses uMain;

procedure SetSize(N, M: Integer; var R: TMatrix);

var

I, J: Integer;

begin

//Установим число строк матрицы (+1 т.к. нумерация идет от 0)

SetLength(R, N+1);

for I := 1 to N do

begin

//Установим число столбцов для каждой строки(+1 т.к. нумерация идет от 0)

SetLength(R[I], M+1);

//Заполним строку нулями

for J := 1 to M do

R[I, J] := 0;

end;

end;

procedure GetSize(R: TMatrix; var N, M: Integer);

begin

//Вернем число строк

N := High(R);

//Вернем число столбцов в первой строке 

M := High(R[1]);

end;

procedure MultMatrix(A, B: TMatrix; var R: TMatrix);

var

aN, aM, bN, bM, I, J, K: Integer;

V: Extended;

begin

//Получим размеры матриц

getSize(A, aN, aM);

getSize(B, bN, bM);

//Проверим согласование внутренних размеров

assert(aM = bN);

//Установим размер результирующей матрицы

SetSize(aN, bM, R);

//Выполним умножение матриц

for I := 1 to aN do

for J := 1 to bM do

begin

V := 0;

for K := 1 to aM do

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

R[I, J] := V;

end;

end;

procedure Transpose(A: TMatrix; var R: TMatrix);

var

N, M, I, J: Integer;

begin

//Получим размер матрицы

GetSize(A, N, M);

//Установим размер результата

SetSize(M, N, R);

//Заполним результат 

for I := 1 to N do

for J := 1 to M do

R[J, I] := A[I, J];

end;

procedure zeidel(N:Integer; e:Extended; A:TMatrix; B:TVector; var X:TVector);

var i,j:longint;

s,v,m:Extended;

cnt: Integer;

Base: Extended;

LastM: Extended;

begin

//Инициализируем счетчик итераций

cnt := 0;

//Инициализируем оценку погрешности

Base := 1;

LastM := 1e6;

// Сам алгоритм

for i:=1 to n do

x[i] := 0;

repeat

m:=0;

for i:=1 to n do

begin

s:=0;

for j:=1 to n do

if i<>j then

s:=s+a[i,j]*x[j];