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
MatCom;
end;
procedure FreeMatrix(var Matrix:TMatrix);
//var i:Integer;
begin// for i:=0 to High(Matrix) do SetLength(Matrix
,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,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,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,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,L);
SetLength(Value,H,H)
// for i:=0 to L do SetLength(Value,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,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,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,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,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;
Iy:=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];
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:=Xi/A[i,i];
end;
for i:=0 to pred(DimMat) do X[Jx]:=PVX;
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:=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;
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-X2) 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:=X[i,0];
for i:=0 to Dim-1 do Mid:=A[i,i];
K:=0;
repeat
for i:=0 to Dim-1 do X1:=X2;
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:=X1+W/Mid*(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;
end else MatrixErrStr:='经过'+IntToStr(N)+'次迭代达不到精度';
SetLength(X1,0)
SetLength(X2,0)
SetLength(Mid,0);
end;
end.
