# 矩阵指数

Moler, Cleve, and Charles Van Loan.“Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later.”SIAM Review 45, no. 1 (January 2003):3–49. https://doi.org/10.1137/S00361445024180.

`A = [0 1 2; 0.5 0 1; 2 1 0]`
```A = 3×3 0 1.0000 2.0000 0.5000 0 1.0000 2.0000 1.0000 0 ```
`Asave = A;`

### 方法 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` 使用矩阵指数的经典定义，表示为幂级数

`${e}^{A}=\sum _{k=0}^{\infty }\frac{1}{k!}{A}^{k}.$`

${\mathit{A}}^{0}$ 是与 $\mathit{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` 假定矩阵包含一组完整的特征向量 $\mathit{V}$，使得 $\mathit{A}=\mathrm{VD}{\mathit{V}}^{-1}$。矩阵指数可以通过对特征值的对角矩阵求幂来计算：

`${e}^{A}=V{e}^{D}{V}^{-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-13 × -0.0711 -0.0444 -0.0799 -0.0622 -0.0488 -0.0933 -0.0711 -0.0533 -0.1066 ```

### 泰勒级数失败

```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 ```

### 特征值和特征向量失败

```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 ```