Computes the exponential of a matrix.
B = linalg_expm(A)
a n-by-n matrix of doubles
a n-by-n matrix of doubles, B=exp(A)
Scales the matrix by power of 2 to increase convergence. Uses Padé approximants.
A = [1 2; 3 4]; expected = [ 5.196895619870500838D+01 7.473656456700319950D+01 1.121048468505048135D+02 1.640738030492098005D+02 ]; F = linalg_expm(A); ulperror = ceil(abs(F-expected)./abs(expected)/%eps) // Test #2 : a 3x3 matrix A = matrix((1:9),3,3); expected = [ 1.118906699413182214D+06 2.533881041898961645D+06 3.948856384384742938D+06 1.374815062935801689D+06 3.113415031380542088D+06 4.852012999825284816D+06 1.630724426458421396D+06 3.692947020862122998D+06 5.755170615265826695D+06 ]; F = linalg_expm(A); ulperror = ceil(abs(F-expected)./abs(expected)/%eps) // Test #3 : from wikipedia A = [ 21 17 6 -5 -1 -6 4 4 16 ]; expected = 1/4 * [ 13*exp(16)-exp(4) 13*exp(16)-5*exp(4) 2*exp(16)-2*exp(4) -9*exp(16)+exp(4) -9*exp(16)+5*exp(4) -2*exp(16)+2*exp(4) 16*exp(16) 16*exp(16) 4*exp(16) ]; F = linalg_expm(A); ulperror = ceil(abs(F-expected)./abs(expected)/%eps) | ![]() | ![]() |
"Matrix Computations", Golub and Van Loan, Third Edition, The Johns Hopkins University Press, p573, Algorithm 11.3-1.
http://en.wikipedia.org/wiki/Matrix_exponential