这个例子展示了计算矩阵指数的19种方法中的3种。有关矩阵指数计算的背景资料，请参阅“25年后计算矩阵指数的19种可疑方法”，SIAM评论45,3-49,2003。我们亦强烈建议使用「伪谱闸道」。本网站的网址为:http://web.comlab.ox.ac.uk/projects/pseudospectra/工具/原料more电脑matlab软件方法/步骤1/14分步阅读

从矩阵A开始

A =

Asave = A;

2/14

按“Enter”键。

得图1所示。

3/14

缩放和平方

是算法11.3.1在Golub和Van Loan中的实现，矩阵计算，第3版。

% Scale A by power of 2 so that its norm is 1/2 .

= 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

4/14

按“Enter”键。

得图2所示。

5/14

exmdem2使用了矩阵指数的经典定义。

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

6/14

按“Enter”键。

得图3所示。

7/14

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

expmdemo3假设矩阵具有完整的特征向量集。作为一种实用的数值方法，精度取决于特征向量矩阵的条件。

A = Asave;

= eig(A);

E = V * diag(exp(diag(D))) / V;

E3 = E

8/14

按“Enter”键。

得图4所示。

9/14

比较结果

对于这个矩阵，它们都做得同样好。

E = expm(Asave);

err1 = E - E1

err2 = E - E2

err3 = E - E3

10/14

按“Enter”键。

得图5所示。

11/14

泰勒级数失败

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

A = ;

E1 = expmdemo1(A)

E2 = expmdemo2(A)

E3 = expmdemo3(A)

12/14

按“Enter”键。

得图6所示。

13/14

这是一个没有完整特征向量集合的矩阵，因此，expmdemo3失败了。

A = ;

E1 = expmdemo1(A)

E2 = expmdemo2(A)

E3 = expmdemo3(A)

14/14

按“Enter”键。

得图7所示。