矩阵指数

打开实时脚本

这个例子展示了计算矩阵指数的19种方法中的3种。

有关矩阵指数计算的背景知识，请参见:

莫勒，克利夫和查尔斯·范洛安。《25年后计算矩阵指数的19种可疑方法》暹罗审查45岁的没有。1(2003年一月):3-49。https://doi.org/10.1137/S00361445024180．

首先创建一个矩阵一个．

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

一个=3×30 1.0000 2.0000 0.5000 0 1.0000 2.0000 1.0000 0

Asave = A;

方法一:缩放和平方

expmdemo1是书中算法11.3.1的实现:

Golub, Gene H.和Charles Van Loan。矩阵计算，第三版。马里兰州巴尔的摩:约翰·霍普金斯大学出版社，1996年。

将A按2的幂进行缩放，使其范数< 1/2。[f,e] = log2(范数(A，“正”));S = max(0,e+1);A = A/2^s;exp的% Pade近似(A)X = a;C = 1/2;E =眼睛(大小(A)) + c*A;D =眼睛(大小(A)) - c*A;Q = 6;P = 1;为K = 2:q c = c *(q- K +1) / (K *(2*q- K +1));X = a *X;cX = c*X;E = E + cX;如果p D = D + cX;其他的D = D - cX;结束P = ~ P;结束E = d \E;通过重复平方撤消缩放为k = 1:s E = E*E;结束E1 = e

E1 =3×35.3091 4.0012 5.5778 2.8088 2.8845 3.1930 5.1737 4.0012 5.7132

方法二:泰勒级数

expmdemo2使用幂级数给出的矩阵指数的经典定义

$e^{一个} ＝ \sum_{k ＝ 0}^{\infty} \frac{1}{k ！} {一个}^{k} ．$

${一个}^{0}$ 单位矩阵的维数和 $一个$ ．作为一种实用的数值方法，该方法速度慢且不准确规范(一)太大了。

A = Asave;exp的%泰勒级数(A)E = 0(大小(A));F =眼睛(大小(A));K = 1;而范数(E+F-E,1) > 0 E = E+F;F = A*F/k;K = K +1;结束E2 = e

E2 =3×35.3091 4.0012 5.5778 2.8088 2.8845 3.1930 5.1737 4.0012 5.7132

方法三:特征值和特征向量

expmdemo3假设矩阵有一个完整的特征向量集合 $V$ 这样 $一个＝ VD V^{- 1}$ ．矩阵指数可以通过对特征值对角线矩阵求幂来计算:

$e^{一个} ＝ V e^{D} V^{- 1} ．$

作为一种实用的数值方法，其精度取决于特征向量矩阵的条件。

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

E3 =3×35.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×310^－１４× 0.3553 0.1776 0.0888 0.0888 0.1332 -0.0444 0 0 -0.2665

err2 = E - E2

err2 =3×310^－１４× 0 0 -0.1776 -0.0444 0 -0.0888 0.1776 0 0.0888

err3 = E - E3

err3 =3×310^－１４× -0.7105 -0.5329 -0.7105 -0.6217 -0.5773 -0.8882 -0.6217 -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×210⁶× -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×20.3679 0.3679 0 0.3679

E2 = expmdemo2(A)

E2 =2×20.3679 0.3679 0 0.3679

E3 = expmdemo3(A)

E3 =2×20.3679 0 0 0.3679

另请参阅

expm