Spencer Graves wrote:I think you need the Schur decomposition, which seems currently not to be available in R. The documentation for the the Matrix package describes a "Schur" function, but it's not available in the current Matrix package, as mentioned in my post on this issue yesterday (subj: Schur decomposition).
Moreover, Lindsey's mexp in rmutils won't work either:
> A = matrix(cbind(c(-1,1),c(-4,3)),nrow=2)
> mexp(A)
Error in solve.default(z$vectors) : system is computationally singular: reciprocal condition number = 4.13756e-017
> mexp(A, "series")
Error in t * x : non-numeric argument to binary operator
>
Have you considered adding a little noise: > mexp(A+1e-6*rnorm(4))
[,1] [,2]
[1,] -2.718284 -10.873139
[2,] 2.718284 8.154859
> mexp(A+1e-6*rnorm(4))
[,1] [,2]
[1,] -2.718284-1.060041e-12i -10.873126-1.575184e-12i
[2,] 2.718284+7.492895e-13i 8.154846+1.578515e-12i
> mexp(A+1e-6*rnorm(4))
[,1] [,2]
[1,] -2.718284+6.146195e-13i -10.873127+1.586731e-12i
[2,] 2.718284-4.004574e-13i 8.154847-8.515411e-13i
> mexp(A+1e-6*rnorm(4))
[,1] [,2]
[1,] -2.718283+3.077538e-13i -10.873130+2.433609e-13i
[2,] 2.718283-3.197442e-14i 8.154847-2.782219e-13i
>
hope this helps. spencer graves
(I believe this is described in one of Richard Bellman's matrix analysis book.)
I rather suspected that that original question was about the matrix exponential and linear systems of differential equations. The solution to such a system can only be written using the matrix exponential for diagonalizable systems and this is the classic example of a nondiagonalizable system.
Calculation of the matrix exponential is an operation that seems straightforward in theory and can be very difficult in practice. Moler and van Loan have a classic paper on "Nineteen Dubious Ways to Calculate the Matrix Exponential" which I would recommend reading.
interesting, I would not have suspected such difficulties based on the technical definition of the matrix exponential I know
exp(A) = Identity + A + A^2/2! + A^3/3! ...
which does not at all require A to be diagonalizable.
Now, the original question was indeed how to solve numerically a linear ODE system with non-diagonalizable matrix. The example I gave in another reply suggests that this could be done by taking all linearly independent eigenvectors, complement them to a full base P, then compute D=invP %*% A %*% P (upper triangular matrix), decompose it into a diagonal matrix + an upper triangular matrix with 0 on the diagonal for which exp(D) would be rather straightforward to compute as long as the system is not too large.
What I do not know is whether the complementing to a full base is feasible. Any idea?
Thanks, Christian.
-- *********************************************************** http://cognition.ups-tlse.fr/vas-y.php?id=chj [EMAIL PROTECTED] Christian Jost (PhD, MdC) Centre de Recherches sur la Cognition Animale Universite Paul Sabatier, Bat IV R3 118 route de Narbonne 31062 Toulouse cedex 4, France Tel: +33 5 61 55 64 37 Fax: +33 5 61 55 61 54
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
