矩阵的计算方法
孙悟空故事睡前故事-移印
2023年2月21日发(作者：绝对值教学设计)对称矩阵的计算方法
1、通过传统方法，把矩阵[AI]通过行变换得到[IA-1]，具体方法其实就是高斯消元法。
2、对矩阵进行LDL'分解，于是A-1=L'-1D-1L-1
其中L-1就是把L的对角线下元素取负即可，D-1就是对角线元素取倒数，L'-1=(L-1)'
第2种方法对对称阵的效果会更好一些。
下面给一段程序（还没有优化）：
type
TMatrix=arrayofarrayofExtended;
ENotSquare=class(Exception);
EDimNotMatch=class(Exception);
...
procedureLDLT(m:TMatrix;varl,d:TMatrix);//LDL'分解
var
n,i,j,k:Integer;
s:Extended;
begin
n:=High(m);
ifHigh(m[0])nthen
('LDL''分解需要方阵！');
SetLength(l,n+1,n+1);
SetLength(d,n+1,n+1);
fori:=0tondo
forj:=0tondobegin
l[i,j]:=Ord(i=j);
d[i,j]:=0;
end;
fori:=0tondobegin
s:=m[i,i];
forj:=0toi-1do
s:=s-Sqr(l[i,j])*d[j,j];
d[i,i]:=s;
fork:=i+1tondobegin
s:=m[k,i];
forj:=0toi-1do
s:=s-l[k,j]*l[i,j]*d[j,j];
l[k,i]:=s/d[i,i];
end;
end;
end;
procedureTranspose(m:TMatrix;varTrans:TMatrix);//矩阵转置
var
i,j,n0,n1:Integer;
begin
n0:=High(m);
n1:=High(m[0]);
SetLength(Trans,n1+1,n0+1);
fori:=0ton0do
forj:=0ton1do
Trans[j,i]:=m[i,j];
end;
procedureMultiply(m1,m2:TMatrix;varRes:TMatrix);//矩阵乘法
var
i,j,k,n0,n1,n2:Integer;
s:Extended;
begin
n0:=High(m1);
n1:=High(m1[0]);
ifn1High(m2)then
('矩阵维数不匹配！');
n2:=High(m2[0]);
SetLength(Res,n0+1,n2+1);
fori:=0ton0do
forj:=0ton2dobegin
s:=0;
fork:=0ton1do
s:=s+m1[i,k]*m2[k,j];
Res[i,j]:=s;
end;
end;
procedureInverseSymmetry(m:TMatrix;varInv:TMatrix);//对称阵求逆
var
l,d,lt,temp:TMatrix;
i,j,k,n:Integer;
s:Extended;
begin
n:=High(m);
LDLT(m,l,d);
SetLength(lt,n+1,n+1);
fori:=0tondo
forj:=0tondo
lt[i,j]:=Ord(i=j);
fori:=0tondobegin
d[i,i]:=1/d[i,i];
forj:=i+1tondobegin
s:=0;
fork:=0toj-1do
s:=s+l[j,k]*lt[k,i];
lt[j,i]:=-s;
end;
end;
Transpose(lt,l);
Multiply(l,d,temp);
Multiply(temp,lt,Inv);
end;