Решение СЛАУ, вычисление обратных матриц и определителей с использованием LU-разложения

Previous  Top  Next

    
 

 

Code:

unit LU_Utils;

 

interface

uses Classes;

type

TDblArr = array of Double;

TIntArr = array of Integer;

TDblMatr = array of array of Double;

 

procedure CopyMatr(Src, Dest: TDblMatr);

 

//исходная матрица разрушается, при необходимости нужно сделать копию

procedure LUDecomp(Matr: TDblMatr; Indx: TIntArr; var DSign: Integer);

procedure LUBackSubst(Matr: TDblMatr; Indx: TIntArr; Reslt: TDblArr);

 

//решение системы линейных уравнений. XY на входе - правые части,

//на выходе - вектор решений

//LU-разложение матрицы может использоваться многократно с

//разными правыми частями (при модификации кода)

procedure SLAUSolver(Matr: TDblMatr; XY: TDblArr);

 

//исходная матрица заменяется обратной

procedure MatrInverse(Matr: TDblMatr);

 

//исходная матрица сохраняется

function Determinant(Matr: TDblMatr): Double;

 

implementation

uses SysUtils;

 

procedure CopyMatr(Src, Dest: TDblMatr);

var

n, i: Integer;

begin

n := Length(Src);

for i := 0 to n - 1 do

   Move(Src[i, 0], Dest[i, 0], n * SizeOf(Double));

end;

 

procedure LUDecomp(Matr: TDblMatr; Indx: TIntArr; var DSign: Integer);

const

Tiny = 1.0E-40;

var

n, k, j, imax, i: Integer;

sum, dum, big: Double;

VV: TDblArr;

DumPtr: Pointer;

begin

n := Length(Matr);

SetLength(VV, n);

DSign := 1;

 

for i := 0 to n - 1 do begin

   big := 0;

   for j := 0 to n - 1 do

     if (abs(Matr[i, j]) > big) then

       big := Abs(Matr[i, j]);

   if (big = 0.0) then

     raise Exception.Create('Stopped. Matrix is singular!!!');

   VV[i] := 1 / big;

end;

 

for j := 0 to n - 1 do begin

   for i := 0 to j - 1 do begin

     sum := Matr[i, j];

     for k := 0 to i - 1 do

       sum := sum - Matr[i, k] * Matr[k, j];

     Matr[i, j] := sum;

   end;

 

   big := 0;

   for i := j to n - 1 do begin

     sum := Matr[i, j];

     for k := 0 to j - 1 do

       sum := sum - Matr[i, k] * Matr[k, j];

     Matr[i, j] := sum;

     dum := VV[i] * abs(sum);

     if (dum > big) then begin

       big := dum;

       imax := i;

     end;

   end;

 

   if (j <> imax) then begin

     DumPtr := Pointer(Matr[imax]);

     Pointer(Matr[imax]) := Pointer(Matr[j]);

     Pointer(Matr[j]) := DumPtr;

     DSign := -DSign;

     VV[imax] := VV[j]

   end;

 

   Indx[j] := imax;

   if (Matr[j, j] = 0) then

     Matr[j, j] := Tiny;

   if (j <> n - 1) then begin

     dum := 1 / Matr[j, j];

     for i := j + 1 to n - 1 do

       Matr[i, j] := Matr[i, j] * dum;

   end;

end;

end;

 

procedure LUBackSubst(Matr: TDblMatr; Indx: TIntArr; Reslt: TDblArr);

var

n, j, ip, ii, i: Integer;

sum: Double;

begin

n := Length(Matr);

ii := -1;

 

for i := 0 to n - 1 do begin

   ip := Indx[i];

   sum := Reslt[ip];

   Reslt[ip] := Reslt[i];

   if (ii <> -1) then

     for j := ii to i - 1 do

       sum := sum - Matr[i, j] * Reslt[j]

   else

     if (sum <> 0.0) then

       ii := i;

   Reslt[i] := sum;

end;

 

for i := n - 1 downto 0 do begin

   sum := Reslt[i];

   if (i < n - 1) then

     for j := i + 1 to n - 1 do

       sum := sum - Matr[i, j] * Reslt[j];

   Reslt[i] := sum / Matr[i, i];

end

end;

 

procedure SLAUSolver(Matr: TDblMatr; XY: TDblArr);

var

Indx: TIntArr;

DSign: Integer;

begin

SetLength(Indx, Length(XY));

LUDecomp(Matr, Indx, DSign);

LUBackSubst(Matr, Indx, XY);

end;

 

procedure MatrInverse(Matr: TDblMatr);

var

n, i, j: Integer;

Indx: TIntArr;

Col: TDblArr;

DSign: Integer;

Temp: TDblMatr;

begin

n := Length(Matr);

SetLength(Indx, n);

SetLength(Col, n);

SetLength(Temp, n, n);

LUDecomp(Matr, Indx, DSign);

for j := 0 to n - 1 do begin

   for i := 0 to n - 1 do

     Col[i] := 0;

   Col[j] := 1;

   LUBackSubst(Matr, Indx, Col);

   for i := 0 to n - 1 do

     Temp[i, j] := Col[i];

end;

CopyMatr(Temp, Matr);

Temp := nil;

end;

 

function Determinant(Matr: TDblMatr): Double;

var

n, j: Integer;

Indx: TIntArr;

DSign: Integer;

Temp: TDblMatr;

begin

n := Length(Matr);

SetLength(Indx, n);

SetLength(Temp, n, n);

CopyMatr(Matr, Temp);

LUDecomp(Temp, Indx, DSign);

Result := DSign;

for j := 0 to n - 1 do

   Result := Result * Temp[j, j];

end;

 

end.

 

 

Code:

unit Unit1;

 

interface

 

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, LU_Utils, StdCtrls;

 

type

TForm1 = class(TForm)

   Memo1: TMemo;

   Button1: TButton;

   Button2: TButton;

   Button3: TButton;

   procedure Button1Click(Sender: TObject);

   procedure Button2Click(Sender: TObject);

   procedure Button3Click(Sender: TObject);

private

   { Private declarations }

public

   { Public declarations }

end;

 

var

Form1: TForm1;

 

implementation

 

{$R *.dfm}

 

function ShowArr(Arr: TDblArr): string;

var

i: Integer;

begin

Result := '';

for i := 0 to Length(Arr) - 1 do

   Result := Result + FormatFloat('0.0000', Arr[i]) + '  ';

end;

 

procedure ShowMatr(Matr: TDblMatr; Strings: TStrings);

var

n, i, j: Integer;

s: string;

begin

n := Length(Matr);

for i := 0 to n - 1 do begin

   s := '';

   for j := 0 to n - 1 do

     s := s + FormatFloat('0.0000', Matr[i, j]) + '  ';

   Strings.Add(s);

end;

end;

 

procedure TForm1.Button1Click(Sender: TObject);

var

M:TDblMatr;

B:TDblArr;

begin

SetLength(M,3,3);

SetLength(B,3);

M[0, 0] := 0.5;

M[0, 1] := 2;

M[0, 2] := -1;

M[1, 0] := -15;

M[1, 1] := -0.1;

M[1, 2] := 3.2;

M[2, 0] := 0.7;

M[2, 1] := 4.1;

M[2, 2] := -0.2;

B[0] := 11.4;

B[1] := -196.02;

B[2] := 10.22;

Memo1.Clear;

ShowMatr(M, Memo1.Lines);

Memo1.Lines.Add('');

SLAUSolver(M, B);

ShowMatr(M, Memo1.Lines);

Memo1.Lines.Add('');

Memo1.Lines.Add(ShowArr(B));

end;

 

procedure TForm1.Button2Click(Sender: TObject);

var

M:TDblMatr;

begin

SetLength(M,3,3);

M[0, 0] := 1;

M[0, 1] := 1;

M[0, 2] := 1;

M[1, 0] := 1;

M[1, 1] := 2;

M[1, 2] := 1;

M[2, 0] := 1;

M[2, 1] := 1;

M[2, 2] := 3;

Memo1.Clear;

ShowMatr(M, Memo1.Lines);

Memo1.Lines.Add('');

MatrInverse(M);

ShowMatr(M, Memo1.Lines);

end;

 

procedure TForm1.Button3Click(Sender: TObject);

var

M:TDblMatr;

begin

SetLength(M,3,3);

M[0, 0] := 1;

M[0, 1] := 0;

M[0, 2] := 1;

M[1, 0] := 0;

M[1, 1] := 2;

M[1, 2] := 0;

M[2, 0] := 1;

M[2, 1] := 0;

M[2, 2] := 3;

Memo1.Clear;

ShowMatr(M, Memo1.Lines);

Memo1.Lines.Add('');

Memo1.Lines.Add(FormatFloat('0.0000', Determinant(M)));

end;

 

end.

 

Автор: MBo

©Drkb::04251