Matrix Exponential

Computes the exponential of a matrix.

expm(A, np = 128)

logm(A)

Arguments

A

numeric square matrix.

np

number of points to use on the unit circle.

Details

For an analytic function \(f\) and a matrix \(A\) the expression \(f(A)\) can be computed by the Cauchy integral $$f(A) = (2 \pi i)^{-1} \int_G (zI-A)^{-1} f(z) dz$$ where \(G\) is a closed contour around the eigenvalues of \(A\).

Here this is achieved by taking G to be a circle and approximating the integral by the trapezoid rule.

logm is a fake at the moment as it computes the matrix logarithm through taking the logarithm of its eigenvalues; will be replaced by an approach using Pade interpolation.

Another more accurate and more reliable approach for computing these functions can be found in the R package `expm'.

Value

Matrix of the same size as A.

References

Moler, C., and Ch. Van Loan (2003). Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. SIAM Review, Vol. 1, No. 1, pp. 1–46.

N. J. Higham (2008). Matrix Functions: Theory and Computation. SIAM Society for Industrial and Applied Mathematics.

Author

Idea and Matlab code for a cubic root by Nick Trefethen in his “10 digits 1 page” project.

Note

This approach could be used for other analytic functions, but a point to consider is which branch to take (e.g., for the logm function).

See also

expm::expm

Examples

##  The Ward test cases described in the help for expm::expm agree up to
##  10 digits with the values here and with results from Matlab's expm !
A <- matrix(c(-49, -64, 24, 31), 2, 2)
expm(A)
#>            [,1]      [,2]
#> [1,] -0.7357588 0.5518191
#> [2,] -1.4715176 1.1036382
# -0.7357588 0.5518191
# -1.4715176 1.1036382

A1 <- matrix(c(10,  7,  8,  7,
                7,  5,  6,  5,
                8,  6, 10,  9,
                7,  5,  9, 10), nrow = 4, ncol = 4, byrow = TRUE)
expm(logm(A1))
#>      [,1] [,2] [,3] [,4]
#> [1,]   10    7    8    7
#> [2,]    7    5    6    5
#> [3,]    8    6   10    9
#> [4,]    7    5    9   10
logm(expm(A1))
#>          [,1]     [,2]      [,3]      [,4]
#> [1,] 9.999101 7.074248  7.844689  7.111275
#> [2,] 7.074248 4.783791  6.225884  4.843151
#> [3,] 7.844689 6.225884 10.043020  8.947118
#> [4,] 7.111275 4.843151  8.947118 10.057621

##  System of linear differential equations: y' = M y  (y = c(y1, y2, y3))
M <- matrix(c(2,-1,1, 0,3,-1, 2,1,3), 3, 3, byrow=TRUE)
M
#>      [,1] [,2] [,3]
#> [1,]    2   -1    1
#> [2,]    0    3   -1
#> [3,]    2    1    3
C1 <- 0.5; C2 <- 1.0; C3 <- 1.5
t  <- 2.0; Mt <- expm(t * M)
yt <- Mt