Description Usage Arguments Details Value Note Author(s) References See Also Examples

This function computes the exponential of a square matrix
*A*, defined as the sum from *r=0* to infinity of
*A^r/r!*.
Several methods are provided. The Taylor series and Padé
approximation are very importantly combined with scaling and squaring.

1 2 3 4 5 6 |

`x` |
a square matrix. |

`method` |
The versions with "*O" call the
original Fortran code, whereas the first ones call the BLAS-using
and partly simplified newer code. |

`order` |
an integer, the order of approximation to be used, for
the "Pade" and "Taylor" methods. The best value for this depends on
machine precision (and slightly on |

`trySym` |
logical indicating if |

`tol` |
a given tolerance used to check if |

`do.sparseMsg` |
logical allowing a message about sparse to dense
coercion; setting it |

`preconditioning` |
a string specifying which implementation of
Ward(1977) should be used when |

The exponential of a matrix is defined as the infinite Taylor series

*
exp(M) = I + M + M^2/2! + M^3/3! + …*

For the "Pade" and "Taylor" methods, there is an `"accuracy"`

attribute of the result. It is an upper bound for the L2 norm of the
Cauchy error `expm(x, *, order + 10) - expm(x, *, order)`

.

Currently, only algorithms which are *“ R-code only”* accept

`sparseMatrix`

class in package
Matrix), i.e., currently only `"R_Eigen"`

and
`"Higham08"`

.
The matrix exponential of `x`

.

For a good general discussion of the matrix exponential problem, see Moler and van Loan (2003).

The `"Ward77"`

method:

Vincent Goulet [email protected], and Christophe
Dutang, based on code translated by Doug Bates and Martin Maechler
from the implementation of the corresponding Octave function
contributed by A. Scottedward Hodel [email protected].

The `"PadeRBS"`

method:

Roger B. Sidje, see the EXPOKIT reference.

The `"PadeO"`

and `"TaylorO"`

methods:

Marina Shapira (U Oxford, UK) and David Firth (U Warwick, UK);

The `"Pade"`

and `"Taylor"`

methods are slight
modifications of the "*O" ([O]riginal versions) methods,
by Martin Maechler, using BLAS and LINPACK where possible.

The `"hybrid_Eigen_Ward"`

method by Christophe Dutang is a C
translation of `"R_Eigen"`

method by Martin Maechler.

The `"Higham08"`

and `"Higham08.b"`

(current default) were
written by Michael Stadelmann, see `expm.Higham08`

.

The `"AlMohy-Hi09"`

implementation (**R** code interfacing to
stand-alone C) was provided and donated by Drew Schmidt, U. Tennesse.

Ward, R. C. (1977). Numerical computation
of the matrix exponential with accuracy estimate.
*SIAM J. Num. Anal.* **14**, 600–610.

Roger B. Sidje (1998).
EXPOKIT: Software package for computing matrix exponentials.
ACM - Transactions on Mathematical Software **24**(1), 130–156.

Moler, C and van Loan, C (2003). Nineteen dubious ways to compute
the exponential of a matrix, twenty-five years later.
*SIAM Review* **45**, 3–49. At
http://epubs.siam.org/doi/pdf/10.1137/S00361445024180

Awad H. Al-Mohy and Nicholas J. Higham (2009)
A New Scaling and Squaring Algorithm for the Matrix Exponential.
*SIAM. J. Matrix Anal. & Appl.*, **31**(3), 970–989.

The package vignette for details on the algorithms and calling the function from external packages.

`expm.Higham08`

for `"Higham08"`

.

`expAtv(A,v,t)`

computes *e^{At} v* (for scalar
*t* and *n*-vector *v*) *directly* and more
efficiently than computing *e^{At}*.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 | ```
x <- matrix(c(-49, -64, 24, 31), 2, 2)
expm(x)
expm(x, method = "AlMohy-Hi09")
## ----------------------------
## Test case 1 from Ward (1977)
## ----------------------------
test1 <- t(matrix(c(
4, 2, 0,
1, 4, 1,
1, 1, 4), 3, 3))
expm(test1, method="Pade")
## Results on Power Mac G3 under Mac OS 10.2.8
## [,1] [,2] [,3]
## [1,] 147.86662244637000 183.76513864636857 71.79703239999643
## [2,] 127.78108552318250 183.76513864636877 91.88256932318409
## [3,] 127.78108552318204 163.67960172318047 111.96810624637124
## -- these agree with ward (1977, p608)
## Compare with the naive "R_Eigen" method:
try(
expm(test1, method="R_Eigen")
) ## platform depently, sometimes gives an error from solve
## or is accurate or one older result was
## [,1] [,2] [,3]
##[1,] 147.86662244637003 88.500223574029647 103.39983337000028
##[2,] 127.78108552318220 117.345806155250600 90.70416537273444
##[3,] 127.78108552318226 90.384173332156763 117.66579819582827
## -- hopelessly inaccurate in all but the first column.
##
## ----------------------------
## Test case 2 from Ward (1977)
## ----------------------------
test2 <- t(matrix(c(
29.87942128909879, .7815750847907159, -2.289519314033932,
.7815750847907159, 25.72656945571064, 8.680737820540137,
-2.289519314033932, 8.680737820540137, 34.39400925519054),
3, 3))
expm(test2, method="Pade")
## [,1] [,2] [,3]
##[1,] 5496313853692357 -18231880972009844 -30475770808580828
##[2,] -18231880972009852 60605228702227024 101291842930256144
##[3,] -30475770808580840 101291842930256144 169294411240859072
## -- which agrees with Ward (1977) to 13 significant figures
expm(test2, method="R_Eigen")
## [,1] [,2] [,3]
##[1,] 5496313853692405 -18231880972009100 -30475770808580196
##[2,] -18231880972009160 60605228702221760 101291842930249376
##[3,] -30475770808580244 101291842930249200 169294411240850880
## -- in this case a very similar degree of accuracy.
##
## ----------------------------
## Test case 3 from Ward (1977)
## ----------------------------
test3 <- t(matrix(c(
-131, 19, 18,
-390, 56, 54,
-387, 57, 52), 3, 3))
expm(test3, method="Pade")
## [,1] [,2] [,3]
##[1,] -1.5096441587713636 0.36787943910439874 0.13533528117301735
##[2,] -5.6325707997970271 1.47151775847745725 0.40600584351567010
##[3,] -4.9349383260294299 1.10363831731417195 0.54134112675653534
## -- agrees to 10dp with Ward (1977), p608.
expm(test3, method="R_Eigen")
## [,1] [,2] [,3]
##[1,] -1.509644158796182 0.3678794391103086 0.13533528117547022
##[2,] -5.632570799902948 1.4715177585023838 0.40600584352641989
##[3,] -4.934938326098410 1.1036383173309319 0.54134112676302582
## -- in this case, a similar level of agreement with Ward (1977).
##
## ----------------------------
## Test case 4 from Ward (1977)
## ----------------------------
test4 <-
structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-10,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0),
.Dim = c(10, 10))
attributes(expm(test4, method="Pade"))
max(abs(expm(test4, method="Pade") - expm(test4, method="R_Eigen")))
##[1] 8.746826694186494e-08
## -- here mexp2 is accurate only to 7 d.p., whereas mexp
## is correct to at least 14 d.p.
##
## Note that these results are achieved with the default
## settings order=8, method="Pade" -- accuracy could
## presumably be improved still further by some tuning
## of these settings.
##
## example of computationally singular matrix
##
m <- matrix(c(0,1,0,0), 2,2)
try(
expm(m, m="R_Eigen")
)
## error since m is computationally singular
expm(m, m="hybrid")
## hybrid use the Ward77 method
``` |

```
Loading required package: Matrix
Attaching package: 'expm'
The following object is masked from 'package:Matrix':
expm
[,1] [,2]
[1,] -0.7357588 0.5518191
[2,] -1.4715176 1.1036382
[,1] [,2]
[1,] -0.7357588 0.5518191
[2,] -1.4715176 1.1036382
[,1] [,2] [,3]
[1,] 147.8666 183.7651 71.79703
[2,] 127.7811 183.7651 91.88257
[3,] 127.7811 163.6796 111.96811
attr(,"accuracy")
[1] 9.942847e-12
[,1] [,2] [,3]
[1,] 147.8666 183.7651 71.79703
[2,] 127.7811 183.7651 91.88257
[3,] 127.7811 163.6796 111.96811
[,1] [,2] [,3]
[1,] 5.496314e+15 -1.823188e+16 -3.047577e+16
[2,] -1.823188e+16 6.060523e+16 1.012918e+17
[3,] -3.047577e+16 1.012918e+17 1.692944e+17
attr(,"accuracy")
[1] 53010
[,1] [,2] [,3]
[1,] 5.496314e+15 -1.823188e+16 -3.047577e+16
[2,] -1.823188e+16 6.060523e+16 1.012918e+17
[3,] -3.047577e+16 1.012918e+17 1.692944e+17
[,1] [,2] [,3]
[1,] -1.509644 0.3678794 0.1353353
[2,] -5.632571 1.4715178 0.4060058
[3,] -4.934938 1.1036383 0.5413411
attr(,"accuracy")
[1] 4.375779e-11
[,1] [,2] [,3]
[1,] -1.509644 0.3678794 0.1353353
[2,] -5.632571 1.4715178 0.4060058
[3,] -4.934938 1.1036383 0.5413411
$dim
[1] 10 10
$accuracy
[1] 1.440257e-16
[1] 1.391826e-08
Error in solve.default(V) :
system is computationally singular: reciprocal condition number = 1.00208e-292
[,1] [,2]
[1,] 1 0
[2,] 1 1
```

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.