/usr/lib/R/site-library/expm/demo/exact-fn.R is in r-cran-expm 0.99-1.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 107 108 109 110 111 112 113 114 115 116 117 118 119 | #### "demo": Make these (function) definitions easily available to useRs
#### ---- --> we use them in tests in ../tests/exact-ex.R
## ~~~~~~~~~~~~~~~~~~~
### For nilpotent matrices A, exp(A) is polynomial in A
### Mathworld gives the example of the general 3 x 3 upper triangle
nilA3 <- function(x,y,z) {
## Purpose: simple nilpotent matrix 3x3 A (with A^n = 0 for n >= 3)
## / 0 x z \
## A = [ 0 0 y ]
## \ 0 0 0 /
## and its exact matrix exponential
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 19 Jan 2008
stopifnot((n <- length(x)) == 1, length(y) == 1, length(z) == 1,
is.numeric(x), is.numeric(y), is.numeric(z))
list(A = cbind(0, rbind(matrix(c(x,0,z,y), 2,2), 0)),
expA = cbind(c(1,0,0), c(x,1,0), c(z + x*y/2, y, 1)))
}
## Note that relErr() is copy-pasted from ../inst/test-tools.R :
##
## The relative error typically returned by all.equal:
relErr <- function(target, current)
mean(abs(target - current)) / mean(abs(target))
facMat <- function(n, R_FUN, ev = R_FUN(n), M = rMat(n, R_FUN = R_FUN))
{
## Purpose: Construct random matrix x of which we "know" expm(x)
## because we set x := M %*% diag(ev) %*% solve(M)
## ----------------------------------------------------------------------
## Arguments: n: dimension of matrices
## R_FUN: random number generator function (n)
## ev: numeric length-n vector of eigenvalues
## M: n x n matrix. Note that the default,
## rMat() will give matrices ``not close to singular''
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: Feb 2008
R_FUN <- match.fun(R_FUN)
stopifnot(n > 0, is.numeric(ev), length(ev) == n,
dim(M) == c(n,n), is.numeric(M))
iM <- solve(M)
## D <- diag(ev); A = M %*% D %*% iM
list(A = M %*% (ev * iM),
expA = M %*% (exp(ev) * iM))
}
### --- The 2x2 example with bad condition , see A3 in ./ex2.R
m2ex3 <- function(eps = 0) {
stopifnot(is.numeric(eps), length(eps) == 1)
A <- rbind(c(-1, 1),
c(eps^2, -1))
I.e <- 1 - eps^2 / 2
V <- I.e* rbind( c(-1, 1),
eps*c( 1, 1))
D <- c(-1-eps, -1+eps)
iV <- ## solve(V)
rbind(c(-1, 1/eps),
c( 1, 1/eps)) / (2 * I.e)
## NOTE: kappa(V) = condition_number(V) == 1/eps exactly
useTol <- 2e-16 / eps
stopifnot(all.equal(diag(2), V %*% iV, tol=useTol),
all.equal(A, V %*% diag(D) %*% iV, tol=useTol) )
ch.e <- cosh(eps)
sh.e <- sinh(eps)
list(A = A,
expA = exp(-1) *
rbind(c( ch.e, sh.e/eps),
c(sh.e*eps, ch.e )))
}
###---
rnilMat <- function(n, R_FUN = function(n) rpois(n, lambda=5))
{
## random upper triangular (zero-diagonal) nilpotent n x n matrix
m <- matrix(0, n,n)
ut <- upper.tri(m)
R_FUN <- match.fun(R_FUN)
m[ut] <- R_FUN(sum(ut))
m
}
set.seed(17)
m <- rnilMat(10)
if(FALSE)
Matrix(m)
## 10 x 10 sparse Matrix of class "dtCMatrix"
##
## . 3 10 7 3 4 9 5 9 6
## . . 5 4 3 . 5 6 3 6
## . . . 5 7 7 3 7 5 6
## . . . . 3 7 6 8 2 7
## . . . . . 9 5 2 7 6
## . . . . . . 8 5 4 6
## . . . . . . . 5 5 3
## . . . . . . . . 3 5
## . . . . . . . . . 3
## . . . . . . . . . .
## An interesting example, rounded from above {see ../tests/exact-ex.R} :
dN <- 9*7*320 # 20160
EmN <- matrix(c(dN, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3*dN, dN, 0, 0, 0, 0, 0, 0, 0, 0,
352800, 5*dN, dN, 0, 0, 0, 0, 0, 0, 0,
1018080, 332640, 5*dN, dN, 0, 0, 0, 0, 0, 0,
2235240, 786240, 292320, 3*dN, dN, 0, 0, 0, 0, 0,
9368520, 3483480, 1582560, 413280, 181440, dN, 0, 0, 0, 0,
24676176, 9598680, 5073600, 1562400, 826560, 161280, dN, 0,0,0,
43730160, 17451000, 10051440, 3430560, 1955520, 504000,
5*dN, dN, 0, 0,
68438436, 27747480, 16853760, 6036240, 3638880, 1038240,
252000, 3*dN, dN, 0,
119725855, 49165892, 31046760, 11652480, 7198800, 2264640,
614880, 191520, 3*dN, dN),
10, 10)
xct10 <- list(m = m, expm = EmN / dN, expmNum = EmN, expmDen = dN)
|