方法 1:加权平方
expmdemo1 是以下著作中算法 11.3.1 的实现:
Golub, Gene H. and Charles Van Loan.Matrix Computations, 3rd edition.Baltimore, MD:Johns Hopkins University Press, 1996.
% 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 = 3×3
5.3091 4.0012 5.5778
2.8088 2.8845 3.1930
5.1737 4.0012 5.7132
方法 2:泰勒级数
expmdemo2 使用矩阵指数的经典定义,表示为幂级数
eA=∑k=0∞1k!Ak.
A0 是与 A 具有相同维度的单位矩阵。作为一种实用的数值方法,如果 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 = 3×3
5.3091 4.0012 5.5778
2.8088 2.8845 3.1930
5.1737 4.0012 5.7132
方法 3:特征值和特征向量
expmdemo3 假定矩阵包含一组完整的特征向量 V,使得 A=VDV-1。矩阵指数可以通过对特征值的对角矩阵求幂来计算:
eA=VeDV-1.
作为一种实际的数值方法,准确性由特征向量矩阵的条件确定。
A = Asave;
[V,D] = eig(A);
E = V * diag(exp(diag(D))) / V;
E3 = E
E3 = 3×3
5.3091 4.0012 5.5778
2.8088 2.8845 3.1930
5.1737 4.0012 5.7132
比较结果
对于此示例中的矩阵,所有三种方法都同样有效。
E = expm(Asave);
err1 = E - E1
err1 = 3×3
10-14 ×
0.3553 0.1776 0.0888
0.0888 0.1332 -0.0444
0 0 -0.2665
err2 = E - E2
err2 = 3×3
10-14 ×
0 0 -0.1776
-0.0444 0 -0.0888
0.1776 0 0.0888
err3 = E - E3
err3 = 3×3
10-14 ×
-0.7105 -0.5329 -0.7105
-0.6661 -0.5773 -0.8882
-0.7105 -0.7105 -0.9770
泰勒级数失败
对于某些矩阵,泰勒级数中的项在变为零之前变得非常大。因此,expmdemo2 失败。
A = [-147 72; -192 93];
E1 = expmdemo1(A)
E1 = 2×2
-0.0996 0.0747
-0.1991 0.1494
E2 = expmdemo2(A)
E2 = 2×2
106 ×
-1.1985 -0.5908
-2.7438 -2.0442
E3 = expmdemo3(A)
E3 = 2×2
-0.0996 0.0747
-0.1991 0.1494
特征值和特征向量失败
以下是不包含一组完整的特征向量的矩阵。因此,expmdemo3 失败。
A = [-1 1; 0 -1];
E1 = expmdemo1(A)
E1 = 2×2
0.3679 0.3679
0 0.3679
E2 = expmdemo2(A)
E2 = 2×2
0.3679 0.3679
0 0.3679
E3 = expmdemo3(A)
E3 = 2×2
0.3679 0
0 0.3679