给个思路 :
1 。通过传统方法 ,把矩阵 [A I] 通过行变换得到 [I A-1],具体方法其实就是高斯消元法。
2 。对矩阵进行 LDL' 分解 ,于是 A-1 = L'-1 D-1 L-1
其中 L-1 就是把 L 的对角线下元素取负即可 ,D-1 就是对角线元素取倒数 ,L'-1 = (L-1)'
第 2 种方法对对称阵的效果会更好一些。
下面给一段程序 ( 还没有优化 ):
type
TMatrix = array of array of Extended;
ENotSquare = class(Exception);
EDimNotMatch = class(Exception);
...
procedure LDLT(m: TMatrix; var l, d: TMatrix); //LDL' 分解
var
n, i, j, k: Integer;
s: Extended;
begin
n := High(m);
if High(m[0]) <> n then
Raise ENotSquare.Create('LDL'' 分解需要方阵！ ');
SetLength(l, n + 1, n + 1);
SetLength(d, n + 1, n + 1);
for i := 0 to n do
for j := 0 to n do begin
l[i, j] := Ord(i = j);
d[i, j] := 0;
end;
for i := 0 to n do begin
s := m[i, i];
for j := 0 to i - 1 do
s := s - Sqr(l[i, j]) * d[j, j];
d[i, i] := s;
for k := i + 1 to n do begin
s := m[k, i];
for j := 0 to i - 1 do
s := s - l[k, j] * l[i, j] * d[j, j];
l[k, i] := s / d[i, i];
end;
end;
end;
procedure Transpose(m: TMatrix; var Trans: TMatrix); // 矩阵转置
var
i, j, n0, n1: Integer;
begin
n0 := High(m);
n1 := High(m[0]);
SetLength(Trans, n1 + 1, n0 + 1);
for i := 0 to n0 do
for j := 0 to n1 do
Trans[j, i] := m[i, j];
end;
procedure Multiply(m1, m2: TMatrix; var Res: TMatrix); // 矩阵乘法
var
i, j, k, n0, n1, n2: Integer;
s: Extended;
begin
n0 := High(m1);
n1 := High(m1[0]);
if n1 <> High(m2) then
Raise EDimNotMatch.Create(' 矩阵维数不匹配！ ');
n2 := High(m2[0]);
SetLength(Res, n0 + 1, n2 + 1);
for i := 0 to n0 do
for j := 0 to n2 do begin
s := 0;
for k := 0 to n1 do
s := s + m1[i, k] * m2[k, j];
Res[i, j] := s;
end;
end;
procedure InverseSymmetry(m: TMatrix; var Inv: 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);
for i := 0 to n do
for j := 0 to n do
lt[i, j] := Ord(i = j);
for i := 0 to n do begin
d[i, i] := 1 / d[i, i];
for j := i + 1 to n do begin
s := 0;
for k := 0 to j - 1 do
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;
评论