矩阵指数

有关矩阵指数计算的背景，请参阅 "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later," SIAM Review 45, 3-49, 2003。

此外，也强烈建议您了解“Pseudospectra Gateway”。网址为：

http://web.comlab.ox.ac.uk/projects/pseudospectra/

下面列出了 19 种计算矩阵指数的方法中的三种。

首先构造矩阵 A

A = [0 1 2; 0.5 0 1; 2 1 0]
Asave = A;

A =

         0    1.0000    2.0000
    0.5000         0    1.0000
    2.0000    1.0000         0

缩放计算和平方计算

expmdemo1 是 Golub 和 Van Loan 合著的《Matrix Computations》第三版中算法 11.3.1 的实现。

% Scale A by power of 2 so that its norm is < 1/2 .
[f,e] = log2(norm(A,'inf'));
s = max(0,e+1);
A = A/2^s;

% Pade approximation for exp(A)
X = A;
c = 1/2;
E = eye(size(A)) + c*A;
D = eye(size(A)) - c*A;
q = 6;
p = 1;
for k = 2:q
   c = c * (q-k+1) / (k*(2*q-k+1));
   X = A*X;
   cX = c*X;
   E = E + cX;
   if p
     D = D + cX;
   else
     D = D - cX;
   end
   p = ~p;
end
E = D\E;

% Undo scaling by repeated squaring
for k = 1:s
    E = E*E;
end

E1 = E

E1 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

通过泰勒级数计算矩阵指数

expmdemo2 使用矩阵指数的经典定义。作为一种适用数值方法，如果 norm(A) 太大，则计算速度将很慢而且不准确。

A = Asave;

% Taylor series for exp(A)
E = zeros(size(A));
F = eye(size(A));
k = 1;
while norm(E+F-E,1) > 0
   E = E + F;
   F = A*F/k;
   k = k+1;
end

E2 = E

E2 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

通过特征值和特征向量计算矩阵指数

expmdemo3 假定矩阵包含一组完整的特征向量。作为一种适用的数值方法，准确性由特征向量矩阵的条件确定。

A = Asave;

[V,D] = eig(A);
E = V * diag(exp(diag(D))) / V;

E3 = E

E3 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

比较结果

对于以下矩阵，结果都具有较高的准确性

E = expm(Asave);
err1 = E - E1
err2 = E - E2
err3 = E - E3

err1 =

  1.0e-014 *

    0.3553    0.1776         0
    0.0444    0.0444   -0.0444
   -0.0888   -0.0888   -0.3553


err2 =

  1.0e-014 *

         0         0   -0.2665
   -0.0888   -0.0888   -0.0888
    0.0888   -0.0888         0


err3 =

  1.0e-013 *

   -0.0799   -0.0444   -0.0711
   -0.0666   -0.0577   -0.0933
   -0.0799   -0.0622   -0.1155

泰勒级数失败

下面的矩阵中泰勒级数中各项在变为零之前非常大。因此，expmdemo2 失败。

A = [-147 72; -192 93];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)

E1 =

   -0.0996    0.0747
   -0.1991    0.1494


E2 =

  1.0e+006 *

   -1.1985   -0.5908
   -2.7438   -2.0442


E3 =

   -0.0996    0.0747
   -0.1991    0.1494

特征值和特征向量失败

以下是不包含一组完整的特征向量的矩阵。因此，expmdemo3 失败。

A = [-1 1; 0 -1];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)

E1 =

    0.3679    0.3679
         0    0.3679


E2 =

    0.3679    0.3679
         0    0.3679


E3 =

    0.3679         0
         0    0.3679