Main Content

Exponenciales de una matriz

Este ejemplo muestra tres de las 19 maneras de calcular la exponencial de una matriz.

Para obtener más información sobre el cálculo de las exponenciales de una matriz, consulte:

Moler, Cleve y Charles Van Loan. "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later". SIAM Review 45, n.º 1 (enero 2003): 3–49. https://doi.org/10.1137/S00361445024180.

Comience creando una matriz A.

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;

Método 1: Escalado y ajuste

expmdemo1 es una implementación del algoritmo 11.3.1 del libro:

Golub, Gene H. y Charles Van Loan. Matrix Computations, 3.ª edición. 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

Método 2: Serie de Taylor

expmdemo2 utiliza la definición clásica de la exponencial de una matriz dada por la serie de potencias

eA=k=01k!Ak.

A0 es la matriz de identidad con las mismas dimensiones que A. Como método numérico práctico, este enfoque es lento e inexacto si norm(A) es demasiado grande.

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

Método 3: Valores propios y vectores propios

expmdemo3 asume que la matriz cuenta con un conjunto completo de vectores propios V como A=VDV-1. La exponencial de una matriz se puede calcular exponenciando la matriz diagonal de los valores propios:

eA=VeDV-1.

Como método numérico práctico, la exactitud viene determinada por la condición de la matriz de vectores propios.

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

Comparar resultados

En la matriz de este ejemplo, los tres métodos funcionan igual de bien.

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

Error de la serie de Taylor

En algunas matrices, los términos de la serie de Taylor aumentan considerablemente antes de bajar a cero. Por consiguiente, expmdemo2 genera un error.

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

Error de valores propios y vectores propios

A continuación, se muestra una matriz sin un conjunto completo de vectores propios. Por consiguiente, expmdemo3 genera un error.

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

Consulte también