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

Материал из DRKB


LU_Utils.pas[править | править код]

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.


Пример использования[править | править код]

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.


Author: MBo
ID: 04251