Please try the following code:
def pot(D,k):
return diagonal_matrix([l^k for l in D.diagonal()])
var('t p k')
R.<q> = PolynomialRing(QQ, 'q')
M = matrix([(-q + 1, q, 0, 0, 0),
(0, -1/2*q + 1, 1/6*q, 1/3*q, 0),
(0, 0, -2/3*q + 1, 2/3*q, 0),
(0, 0, 1/4*q, -1/2*q + 1, 1/4*q),
(0, 0, 0, 0, 1)])
poly = M.characteristic_polynomial()
R2 = PolynomialRing(QQ, ['x','q'])
poly2 = R2(poly)
print factor(poly2)
P = matrix(sum((f(x=M)(q=t).kernel().basis() for f,expo in factor(poly2
)),[]))
D = P*M*~P
Mk = ~P * pot(D,k) * P
Mk
the goal is: given a matrix whose entries are (linear) polynomials in one
variable 'q', with rational coefficients, find the symbolic function f_ij
that, given q and k, returns the ij-th entry to the matrix M^k.
This comes from a problem in Markov chains. You can find that here:
https://malabares.cancamusa.net/home/pub/12/
it's in spanish, but I'm sure you can guess the relevant info from the
pictures. Idea: take any absorbing Markov chain:
https://en.wikipedia.org/wiki/Absorbing_Markov_chain
and try to find the function f(n) that gives the probability that the state
changes from state 1 to a terminal state. You need a symbolic power of a
rational or real matrix, which is problematic but can be done:
http://ask.sagemath.org/question/709/an-nvarn
However this time it's a symbolic power of a symbolic matrix. The point I
want to make here is that Sage can factor the characteristic polynomial,
and find the Jordan form of the matrix, but for that it's necesary to do
some hocus pocus changing from the symbolic ring to polynomial ring in two
variables, back and forth. The code is fragile: change the letters for the
variables, change the rings, and it doesn't work. And the questions are:
- is this example too particular to care?
- can Sage do better with a little bit of help?
- are there big issues with symbolic power of matrices, or should I send a
patch with my code that works-but-I-dont-know-the-theory?
- do you fancy a wrapper around a graph with edge labels that works better
for Markov chains?
Thanks for your time
--
You received this message because you are subscribed to the Google Groups
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-devel+unsubscr...@googlegroups.com.
To post to this group, send email to sage-devel@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-devel.
For more options, visit https://groups.google.com/groups/opt_out.
