Last data update: 2014.03.03
|
R: Matrix exponential
Matrix exponential
Description
Matrix exponential af a square matrix computed by the pade
approximation.
Usage
matexp(a, dt=1,order = 8)
Arguments
a |
A square numeric matrix
|
dt |
Integration Time step
|
order |
Pade approximation order
|
Details
This implementation is based on Niels Rode Kristensens work. This package is also highly inspired
by David Firth's R package mexp.
Value
The matrix exponential is returned. The function issues an error if
problems occured in the fortran engine.
Note
For indepth material on matrix exponentials - see Moler and van Loan (2003).
Author(s)
S<c3><b8>ren Klim, Stig B. Mortensen
References
This implementation is based on Niels Rode Kristensens work. This package is also highly inspired
by David Firth's R package mexp.
The examples below are all from David Firth's mexp package but the
accuracy example has been removed as this package does not calculate
the accuracy.
Niels Rode Kristensen, http://www2.imm.dtu.dk/~ctsm/
Examples
##
## The test cases have been taken directly from David Firths MEXP package.
##
##
## ----------------------------
## Test case 1 from Ward (1977)
## ----------------------------
test1 <- t(matrix(c(
4, 2, 0,
1, 4, 1,
1, 1, 4), 3, 3))
matexp(test1)
## 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)
##
## A naive alternative to mexp, using spectral decomposition:
mexp2 <- function(matrix){
z <- eigen(matrix,sym=FALSE)
Re(z$vectors %*% diag(exp(z$values)) %*%
solve(z$vectors))
}
try(
mexp2(test1)
) ## now gives an error from solve !
##
## 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))
matexp(test2)
## [,1] [,2] [,3]
##[1,] 5496313853692357 -18231880972009844 -30475770808580828
##[2,] -18231880972009852 60605228702227024 101291842930256144
##[3,] -30475770808580840 101291842930256144 169294411240859072
## -- which agrees with Ward (1977) to 13 significant figures
mexp2(test2)
## [,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))
matexp(test3)
## [,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.
mexp2(test3)
## [,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).
##
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(PSM)
Loading required package: MASS
Loading required package: numDeriv
Loading required package: deSolve
Attaching package: 'deSolve'
The following object is masked from 'package:graphics':
matplot
Loading required package: ucminf
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/PSM/matexp.Rd_%03d_medium.png", width=480, height=480)
> ### Name: matexp
> ### Title: Matrix exponential
> ### Aliases: matexp
> ### Keywords: math
>
> ### ** Examples
>
> ##
> ## The test cases have been taken directly from David Firths MEXP package.
> ##
> ##
> ## ----------------------------
> ## Test case 1 from Ward (1977)
> ## ----------------------------
> test1 <- t(matrix(c(
+ 4, 2, 0,
+ 1, 4, 1,
+ 1, 1, 4), 3, 3))
> matexp(test1)
[,1] [,2] [,3]
[1,] 147.8666 183.7651 71.79703
[2,] 127.7811 183.7651 91.88257
[3,] 127.7811 163.6796 111.96811
> ## 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)
> ##
> ## A naive alternative to mexp, using spectral decomposition:
> mexp2 <- function(matrix){
+ z <- eigen(matrix,sym=FALSE)
+ Re(z$vectors %*% diag(exp(z$values)) %*%
+ solve(z$vectors))
+ }
> try(
+ mexp2(test1)
+ ) ## now gives an error from solve !
[,1] [,2] [,3]
[1,] 147.8666 183.7651 71.79703
[2,] 127.7811 183.7651 91.88257
[3,] 127.7811 163.6796 111.96811
> ##
> ## 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))
> matexp(test2)
[,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,] 5496313853692357 -18231880972009844 -30475770808580828
> ##[2,] -18231880972009852 60605228702227024 101291842930256144
> ##[3,] -30475770808580840 101291842930256144 169294411240859072
> ## -- which agrees with Ward (1977) to 13 significant figures
> mexp2(test2)
[,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,] 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))
> matexp(test3)
[,1] [,2] [,3]
[1,] -1.509644 0.3678794 0.1353353
[2,] -5.632571 1.4715178 0.4060058
[3,] -4.934938 1.1036383 0.5413411
> ## [,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.
> mexp2(test3)
[,1] [,2] [,3]
[1,] -1.509644 0.3678794 0.1353353
[2,] -5.632571 1.4715178 0.4060058
[3,] -4.934938 1.1036383 0.5413411
> ## [,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).
> ##
>
>
>
>
>
>
> dev.off()
null device
1
>
|