注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

鑫淼梦园的博客

圆你的梦想 从这里开始

 
 
 

日志

 
 

对称逆矩阵/矩阵乘的计算方法  

2013-01-27 13:22:06|  分类: delphi xe3 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
给个思路 :
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;

  评论这张
 
阅读(283)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017