Метод 1: Масштабирование и обработка на квадрат

expmdemo1 реализация алгоритма 11.3.1 в книге:

Golub, Джин Х. и Чарльз Ван Лоун. Матричные Расчеты, 3-й выпуск. Балтимор, 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 использует классическое определение для матричного экспоненциала, данного степенным рядом

${\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}$. Матричный экспоненциал может быть вычислен возведением в степень диагональная матрица собственных значений:

Как практический численный метод, точность определяется условием матрицы собственного вектора.

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

Отказ ряда Тейлора

Для некоторых матриц условия в Ряду Тейлора становятся очень большими, прежде чем они перейдут к нулю. Следовательно, 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