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
es la matriz de identidad con las mismas dimensiones que . 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 como . La exponencial de una matriz se puede calcular exponenciando la matriz diagonal de los valores propios:
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-14 ×
-0.7105 -0.5329 -0.7105
-0.6217 -0.5773 -0.8882
-0.6217 -0.7105 -0.9770
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