SCM

SCM Repository

[matrix] Annotation of /pkg/tests/matr-exp.R
ViewVC logotype

Annotation of /pkg/tests/matr-exp.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 532 - (view) (download)

1 : maechler 522 library(Matrix)
2 :    
3 : maechler 519 ## Matrix Exponential
4 :    
5 :     assert.EQ.mat <- function(M, m, tol = 1e-15)
6 :     stopifnot(all.equal(as(M, "matrix"), m, tol = tol))
7 :    
8 :    
9 :     m1 <- Matrix(c(1,0,1,1), nc = 2)
10 :     e1 <- expm(m1)
11 : maechler 521 assert.EQ.mat(e1, cbind(c(exp(1),0), exp(1)))
12 : maechler 519
13 :     m2 <- Matrix(c(-49, -64, 24, 31), nc = 2)
14 :     e2 <- expm(m2)
15 : maechler 532 ## The true matrix exponential is 'te2':
16 :     e_1 <- exp(-1)
17 :     e_17 <- exp(-17)
18 :     te2 <- rbind(c(3*e_17 - 2*e_1, -3/2*e_17 + 3/2*e_1),
19 :     c(4*e_17 - 4*e_1, -2 *e_17 + 3 *e_1))
20 :     assert.EQ.mat(e2, te2, tol = 1e-13)
21 :     ## See the (average relative) difference:
22 :     all.equal(as(e2,"matrix"), te2, tol = 0) # 1.48e-14 on "lynne"
23 :    
24 : maechler 520 ## The ``surprising identity'' det(exp(A)) == exp( tr(A) )
25 :     ## or log det(exp(A)) == tr(A) :
26 :     stopifnot(all.equal(determinant(e2)$modulus, sum(diag(m2))))
27 : maechler 519
28 :     m3 <- Matrix(cbind(0,rbind(6*diag(3),0)), nc = 4)
29 :     e3 <- expm(m3)
30 :     assert.EQ.mat(e3,
31 :     rbind(c(1,6,18,36),
32 :     c(0,1, 6,18),
33 :     c(0,0, 1, 6),
34 :     c(0,0, 0, 1)))
35 :    

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge