方法 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