Matrix.pas in mathematicssmallhelp


tool commonly used mathematical calculations...Original Link
    Sponsored links

			
			unit Matrix;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  ExtCtrls,stdctrls, buttons,  grids, Menus, comctrls;

type
  TMatrix = array of array of Extended;
  TDataToFloatFun=function (Str:String; var Value:Extended):Boolean of Object;

//------------------------------------
Type
 TMatrixS=Array of TMatrix;
var
 MatrixErrStr:String;

procedure FreeMatrix(var Matrix:TMatrix);

//转置矩阵
function RollMatrix(Matrix:TMatrix; var Value:TMatrix):Boolean;
//求方阵行列式
function AbsMatrix(Matrix:TMatrix; var Value:Extended):Boolean;
//
function IsMatrixOdd(Matrix:TMatrix):Boolean;
//求方阵的逆
function AthWart(Matrix:TMatrix; var Value:TMatrix):Boolean;
//将Matrix1与Matrix2相加(Anti=False)或相减(Anti=True)
function AddMatrix(Matrix1,Matrix2:TMatrix; var Value:TMatrix; Anti:Boolean=False):Boolean;//加法
//将Matrix1(M*K)和Matrix2(K*N)相乘得Value(M*N)
function MulMatrix(Matrix1,Matrix2:TMatrix; var Value:TMatrix):Boolean;  //乘法

//用选主元的高斯消元法求A*X=B 的X方程解
function SysLin(A,B:TMatrix; var X:TMatrix):Boolean;
//用松弛迭代法解方程
function SorGS(A,B:TMatrix; Var X:TMatrix; W,V:Extended; N:Integer):Boolean;

implementation
uses Math,ZXDDialogs;

type
 PMatCom=^TMatCom;
 TMatCom=Array of Record
          Matrix:TMatrix;
          FunIndex:Byte;
          NexFun:Byte;
          Fun:PMatCom;
         end;
procedure FreeMatrix(var Matrix:TMatrix);
//var i:Integer;
begin// for i:=0 to High(Matrix) do SetLength(Matrix[i],0);
 SetLength(Matrix,0,0);
end;
    //转置     行列式值  逆          加          减          乘
const _Roll=1; _Abs=2;  _athwart=3; _Add=128+1; _Sub=128+2; _Mul=128+3;
      ReserveStrCount=12;
      ReserveStr:Array [0..ReserveStrCount-1] of String=('ROLLMAT','ABSMAT','ATHWART',
                       '+','-','*','(',')','[',']','{','}');

function RollMatrix(Matrix:TMatrix; var Value:TMatrix):Boolean;
var
 W,H,i,j:Integer;
begin
 Result:=False;
 H:=Length(Matrix); if H=0 then exit;
 W:=Length(Matrix[0]); if W=0 then exit;
 SetLength(Value,H,W); //for i:=0 to H-1 do SetLength(Value[i],W);
 for i:=0 to H-1 do for j:=0 to W-1 do Value[j,i]:=Matrix[i,j];
 Result:=True;
end;

procedure GetMatrix(Matrix:TMatrix; N,R:Integer; var Mat:TMatrix);
var
 i,j:Integer;
begin
 SetLength(Mat,R,R);// for i:=0 to R-1 do SetLength(Mat[i],R);
 for i:=0 to N-1 do for j:=0 to R-1 do Mat[j,i]:=Matrix[j+1,i];
 for i:=N+1 to R do for j:=0 to R-1 do Mat[j,i-1]:=Matrix[j+1,i];
end;
function GetMatValue(Matrix:TMatrix; R:Integer):Extended;
var i:Integer; M:TMatrix;
begin
 if R=0 then begin Result:=Matrix[0,0]; exit; end;
 Result:=0;
 for i:=0 to R do  begin
  GetMatrix(Matrix,i,R,M);
  if i mod 2=0
   then Result:=Result+Matrix[0,i]*GetMatValue(M,R-1)
   else Result:=Result-Matrix[0,i]*GetMatValue(M,R-1);
 end;
 SetLength(M,0,0);// for i:=0 to R-1 do SetLength(M[i],0); SetLength(M,0);
end;
function AbsMatrix(Matrix:TMatrix; var Value:Extended):Boolean;
var
 W,H:Integer;
begin
 Result:=False;
 H:=Length(Matrix); if H=0 then exit;
 W:=Length(Matrix[0]); if W<>H then exit;
 Value:=GetMatValue(Matrix,H-1);
 Result:=True;
end;

procedure GetMatrixS(Matrix:TMatrix; R,M,N:Integer; var Mat:TMatrix);
var
 i,j:Integer;
begin
 for i:=0 to N-1 do  for j:=0 to M-1 do Mat[j,i]:=Matrix[j,i];
 for i:=0 to N-1 do  for j:=M+1 to R do Mat[j-1,i]:=Matrix[j,i];
 for i:=N+1 to R do  for j:=0 to M-1 do Mat[j,i-1]:=Matrix[j,i];
 for i:=N+1 to R do  for j:=M+1 to R do Mat[j-1,i-1]:=Matrix[j,i];
end;
function ConCom(Matrix:TMatrix; var Value:TMatrix):Boolean; //求伴随矩阵
var
 W,H,L,i,j:Integer;
 X:Extended;
 M:TMatrix;
begin
 Result:=False;
 H:=Length(Matrix);    if H=0 then exit;
 W:=Length(Matrix[0]); if W<>H then exit;

 L:=H-1; SetLength(M,L,L);// for i:=0 to L-1 do SetLength(M[i],L);
 SetLength(Value,H,H);    // for i:=0 to L do SetLength(Value[i],H);

 for i:=0 to L do  for j:=0 to L do
  begin
   GetMatrixS(Matrix,L,i,j,M);
   AbsMatrix(M,X);
   if (i+j)mod 2=0 then Value[i,j]:=X
                   else Value[i,j]:=-X;
  end;
 SetLength(M,0,0);// FreeMatrix(M);
 Result:=True;
end;
function AthWart(Matrix:TMatrix; var Value:TMatrix):Boolean;//求逆矩阵
var
 X:Extended;
 i,j,R:Integer;
 Temp:TMatrix;
begin
 Result:=False;
 if AbsMatrix(Matrix,X) then begin
  if X=0 then begin MatrixErrStr:='该方阵是奇异的'; exit; end;
  Result:=True;

  R:=High(Matrix);  if R=0 then begin Value[0,0]:=1/X; exit; end;

  SetLength(Value,R+1,R+1); //for i:=0 to R do SetLength(Value[i],R+1);
  ConCom(Matrix,Temp);
  for i:=0 to R do for j:=0 to R do Value[j,i]:=Temp[i,j]/X;
  FreeMatrix(Temp);
 end else MatrixErrStr:='该矩阵不是方阵';
end;

function AddMatrix(Matrix1,Matrix2:TMatrix; var Value:TMatrix; Anti:Boolean=False):Boolean;//加法
var
 H,W,i,j:Integer;
begin
 Result:=False; MatrixErrStr:='加法运算的矩阵不匹配';
 H:=High(Matrix1);    if H<>High(Matrix2) then exit; if H<0 then exit;
 W:=High(Matrix1[0]); if W<>High(Matrix2[0]) then exit; if W<0 then exit;
 SetLength(Value,H+1,W+1);// for i:=0 to H do SetLength(Value[i],W+1);
 if Anti
    then for i:=0 to H do for j:=0 to W do Value[i,j]:=Matrix1[i,j]-Matrix2[i,j]
    else for i:=0 to H do for j:=0 to W do Value[i,j]:=Matrix1[i,j]+Matrix2[i,j];
 Result:=True;
end;
function MulMatrix(Matrix1,Matrix2:TMatrix; var Value:TMatrix):Boolean;  //乘法
var
 H1,WH,W2,
 i,j,k:Integer;
 X:Extended;
begin
 Result:=False; MatrixErrStr:='乘法运算的矩阵不匹配';
 H1:=High(Matrix1); if H1<0 then exit;
 WH:=High(Matrix1[0]); if WH<0 then exit;

 if (H1=0) and (WH=0) then begin
  WH:=High(Matrix2); if WH<0 then exit;
  W2:=High(Matrix2[0]); if W2<0 then exit;
  Result:=True; SetLength(Value,WH+1,W2+1); // for i:=0 to WH do SetLength(Value[i],W2+1);
  X:=Matrix1[0,0];
  for i:=0 to WH do for j:=0 to W2 do Value[i,j]:=Matrix2[i,j]*X;
 end;

 if WH<>High(Matrix2) then exit;
 W2:=High(Matrix2[0]); if W2<0 then exit;
 Result:=True; SetLength(Value,H1+1,W2+1); //for i:=0 to H1 do SetLength(Value[i],W2+1);
 for i:=0 to H1 do for j:=0 to W2 do begin
  X:=0; for K:=0 to WH do X:=X+Matrix1[i,K]*Matrix2[k,j];
  Value[i,j]:=X;
 end;
end;

type RCoVecE=Array of Integer; RCOVec=Array of Extended;

function TrianguleO(A:TMatrix;var AT:TMatrix;var Jx,Iy:RCOVecE;DimMat:Integer):Boolean;
var i, k     : integer;
    im , jm  : integer;
    ATkk     : extended;
    C        : extended;
    Max,Det  : extended;
    SignDet  : integer;
    ErrorMat : Integer;

  PROCEDURE InitJxIy;  //原序号
  var i:integer;
  begin
    for i:=0 to pred(DimMat) do
    begin
      Jx[i]:=i;
      Iy[i]:=i;
    end;
  end;

  PROCEDURE CopiAAT;
  var i,j:integer;
  begin
    for i:=0 to pred(DimMat) do
    for j:=0 to pred(DimMat) do
      AT[i,j]:=A[i,j];
  end;

  PROCEDURE Reduit;     //消元
  var i,j:integer;
  begin
    ATkk:=AT[k,k];
      for i:=k+1 to pred(DimMat) do
        begin
          C   :=AT[i,k]/ATkk;
          for j:=k+1 to pred(DimMat) do
            AT[i,j]:=AT[i,j]-C*AT[k,j];
          for j:=0 to k do
            if j<>k
              then AT[i,j]:=AT[i,j]-C*AT[k,j]
              else AT[i,j]:=-C;
        end;
  end;

  PROCEDURE ChercheMax(var Max:extended;var im,jm:integer); //寻找最大值
  var i,j : integer;
      M   : extended;
  begin
    Max:=abs(AT[k,k]);
    im:=k;jm:=k;
    for i:=k to pred(DimMat) do
      for j:=k to pred(DimMat) do
        begin
          M:=abs(AT[i,j]);
          if M>Max then
            begin
              Max:=M;
              im:=i;jm:=j;
            end;
        end;
  end;

  PROCEDURE PermuteCol;  //交换列
  var Ind : integer;
      i   : integer;
      x   : extended;
  begin
    for i:=0 to pred(DimMat) do
      begin
        x:=AT[i,k];
        AT[i,k]:=AT[i,jm];
        AT[i,jm]:=x;
      end;
    Ind   :=Jx[k];
    Jx[k] :=Jx[jm];
    Jx[jm]:=Ind;
    SignDet:=-SignDet;
  end;

  PROCEDURE PermuteRow; //交换行
  var Ind : integer;
      j   : integer;
      x   : extended;
  begin
    for j:=0 to pred(DimMat) do
      begin
        x:=AT[k,j];
        AT[k,j]:=AT[im,j];
        AT[im,j]:=x;
      end;
    Ind:=Iy[k];
    Iy[k]:=Iy[im];
    Iy[im]:=ind;
    SignDet:=-SignDet;
  end;

begin
  InitJxIy;
  CopiAAT;
  k:=0;
  SignDet:=1;    ErrorMat:=0;
  while (k<=pred(DimMat)-1) and (ErrorMat=0) do
    begin
      ChercheMax(Max,im,jm);
      if Max<>0
        then
          begin
            if k<>jm then PermuteCol;
            if k<>im then PermuteRow;
            Reduit;
          end
        else begin ErrorMat:=1; MatrixErrStr:='最大值为0'; end;
      k:=k+1;
    end;
  Det:=SignDet;
  for i:=0 to pred(DimMat) do
    Det:=Det*AT[i,i];
  if Det=0 then Result:=False
           else Result:=True;
end;
PROCEDURE ResouTrO(A:TMatrix;var X,Y:RCOVec;Jx,Iy:RCOVecE; DimMat:Integer);
var i, j     : integer;
    PVX      : RCOVec;
    Xi       : extended;

begin
  SetLength(PVX,DimMat);
  for i:=pred(DimMat) downto 0 do
    begin
      Xi:=Y[Iy[i]];
      for j:=0 to i-1 do Xi:=Xi+A[i,j]*Y[Iy[j]];
      for j:=i+1 to pred(DimMat) do Xi:=Xi-A[i,j]*PVX[j];
      PVX[i]:=Xi/A[i,i];
    end;
  for i:=0 to pred(DimMat) do X[Jx[i]]:=PVX[i];
  Setlength(PVX,0);
end;
function SysLin(A,B:TMatrix; var X:TMatrix):Boolean;
var AT       : TMatrix;
    Jx,Iy    : RCOVecE;
    X2,B2    : RCOVec;
    i:Integer;
    DimMat:Integer;
begin
 Result:=False; MatrixErrStr:='不匹配的矩阵';
 DimMat:=Length(A);
 if DimMat=0 then exit; if DimMat<>Length(A[0]) then exit;
 if Length(B)<>DimMat then exit; if High(B[0])<>0 then exit;

 SetLength(AT,DimMat,DimMat);
 SetLength(Jx,DimMat);
 SetLength(Iy,DimMat);
 SetLength(X2,DimMat);
 SetLength(B2,DimMat);
 for i:=0 to DimMat-1 do B2[i]:=B[i,0];

 Result:=True;
 if TrianguleO(A,AT,Jx,Iy,DimMat)
  then begin ResouTrO(AT,X2,B2,Jx,Iy,DimMat);
   SetLength(X,DimMat,1);
   for i:=0 to DimMat-1 do X[i,0]:=X2[i];
  end else begin Result:=False; MatrixErrStr:='方程无唯一解'; end;

 SetLength(X2,0);
 SetLength(B2,0);
 SetLength(Iy,0);
 SetLength(Jx,0);
 SetLength(At,0,0);
end;

function IsMatrixOdd(Matrix:TMatrix):Boolean;
var
 AT:TMatrix;
 Jx,Iy:RCoVecE;
 DimMat:Integer;
begin
 Result:=False;
 DimMat:=Length(Matrix); if DimMat=0 then exit;
 if DimMat<>Length(Matrix[0]) then exit;

  SetLength(AT,DimMat,DimMat); SetLength(Jx,DimMat); SetLength(Iy,DimMat);
 Result:=TrianguleO(Matrix,AT,Jx,Iy,DimMat);
  SetLength(Iy,0); SetLength(Jx,0); SetLength(At,0,0);
end;

function SorGS(A,B:TMatrix; Var X:TMatrix; W,V:Extended; N:Integer):Boolean;
var
 X1,X2,
 Mid:Array of Extended;
 i,j,K,Dim:Integer;
 Y:Extended;

 function IsVOK:Boolean;
 var i:Integer;
 begin
  Result:=False;
  for i:=0 to Dim-1 do if V<Abs(X1[i]-X2[i]) then exit;
  Result:=True;
  MatrixErrStr:='';
 end;

begin
 Result:=False; MatrixErrStr:='不匹配的矩阵';
 Dim:=Length(A); if Dim=0 then exit;
 if Dim<>Length(A[0]) then exit;
 if Dim<>Length(B) then exit; if High(B[0])<>0 then exit;

 MatrixErrStr:='方程无解'; V:=Abs(V);
 for i:=0 to Dim-1 do if A[i,i]=0 then exit;
 SetLength(X1,Dim); SetLength(X2,Dim); SetLength(Mid,Dim);
 for i:=0 to Dim-1 do X2[i]:=X[i,0];
 for i:=0 to Dim-1 do Mid[i]:=A[i,i];

 K:=0;
 repeat
  for i:=0 to Dim-1 do X1[i]:=X2[i];
  for i:=0 to Dim-1 do begin
   Y:=0;
   for j:=0 to i-1 do Y:=Y+A[i,j]*X2[j];
   for j:=i to Dim-1 do Y:=Y+A[i,j]*X1[j];
   X2[i]:=X1[i]+W/Mid[i]*(B[i,0]-Y);
  end;
  Inc(K);
 until (K>N)or IsVOK;

 if MatrixErrStr='' then begin
  Result:=True;  for i:=0 to Dim-1 do X[i,0]:=X2[i];
 end  else MatrixErrStr:='经过'+IntToStr(N)+'次迭代达不到精度';

 SetLength(X1,0); SetLength(X2,0); SetLength(Mid,0);
end;


end.
			click here to download the whole source code package.

			
			


Project Files

    Sponsored links
NameSizeDate
 <数学小帮手>0.00 B11-01-05 09:31
 AutoSaveDesk.dcr480.00 B18-09-01 11:19
 AutoSaveDesk.pas14.46 kB18-06-02 21:17
 DataMaker.pas44.69 kB18-06-02 21:25
 HighLightedMemo.dcr488.00 B30-07-01 16:05
 HighlightedMemo.pas8.69 kB18-06-02 11:24
 HSBitMapSD.pas14.07 kB18-06-02 21:25
 Matrix.pas11.63 kB18-06-02 21:26
 PackageZXD.cfg286.00 B18-06-02 11:27
 PackageZXD.dof1.03 kB18-06-02 11:27
 PackageZXD.dpk786.00 B18-06-02 11:26
 PackageZXD.res1.50 kB18-06-02 11:25
 PScMain.cfg386.00 B18-06-02 21:28
 PScMain.dof1.07 kB18-06-02 21:28
 PScMain.dpr202.00 B21-03-02 16:30
 PScMain.res876.00 B18-06-02 21:28
 ScienceDraw.dcr480.00 B30-03-01 10:44
 ScienceDraw.pas48.17 kB18-06-02 11:21
 UnitMain.ddp51.00 B18-06-02 11:27
 UnitMain.dfm189.78 kB18-06-02 11:27
 UnitMain.pas98.32 kB18-06-02 11:27
 ZXDDialogs.pas4.58 kB18-06-02 21:28
...

Related Items

    Sponsored links