# Limit of a matrix power

Hello!

I'm new in Sagemath and I want to calculate the limit of the power of a given 3x3 diagonal matrix:

J=matrix([[1,0,0],[0,-1/2-1/2*I,0],[0,0,-1/2+1/2*I]]); #Given diagonal matrix
n=var('n');
limit(J^n,n=infinity) #The limit command fails =(


Is there any way to do this?

Trank you so much!

edit retag close merge delete

Yes, in the past it was! Now I'm using it as an exercise to learn a bit of Sage

I've edited the question to focus on the limit of the diagonal matrix. In Mathematica there is a way to do this directly, but finally I've found a workaround in Sage:

diagonal_matrix([limit((J^n)[k,k],n=infinity) for k in range(0,3)])


Maybe it's not the most efficient way to do this, but this works with diagonal matrices... The difference is that one needs to execute component by component, also for a more general case (non diagonal matrices).

Sort by » oldest newest most voted

Computing the n-th power of a matrix, where n is a symbolic variable, is possible.

Using that, you could solve your problem as follows.

sage: J = matrix([[1, 0, 0], [0, -1/2-1/2*I, 0], [0, 0, -1/2 + 1/2*I]])
sage: n = SR.var('n')
sage: J^n
[               1                0                0]
[               0 (-1/2*I - 1/2)^n                0]
[               0                0  (1/2*I - 1/2)^n]
sage: J.parent()([lim(x, n=infinity) for x in (J^n).list()])
[1 0 0]
[0 0 0]
[0 0 0]


Edit: this is simpler and more readable:

sage: J = matrix([[1, 0, 0], [0, -1/2-1/2*I, 0], [0, 0, -1/2 + 1/2*I]])
sage: n = SR.var('n')
sage: (J^n).apply_map(lambda x: lim(x, n=oo))
[1 0 0]
[0 0 0]
[0 0 0]


You could then define a function and apply it to any matrix.

sage: def limit_of_power(m):
....:     n = SR.var('n')
....:     return (m^n).apply_map(lambda x: limit(x, n=oo))
....:
sage: limit_of_power(J)
[1 0 0]
[0 0 0]
[0 0 0]

more

Really good solution! That's exactly what I was looking for!

With Sage, you can have a look at the norms of the eigenvalues:

sage: [e.norm() for e in M.eigenvalues()]
[1, 0.50000000000000000?, 0.50000000000000000?]


You can have a look at the (right) eigenvectors:

sage: M.eigenvectors_right()
[(1, [
(1, 1/10, 1/20)
], 1),
(-0.50000000000000000? - 0.50000000000000000?*I,
[(1, -0.10000000000000000? + 0.10000000000000000?*I, 0.?e-20 - 0.10000000000000000?*I)],
1),
(-0.50000000000000000? + 0.50000000000000000?*I,
[(1, -0.10000000000000000? - 0.10000000000000000?*I, 0.?e-20 + 0.10000000000000000?*I)],
1)]


This should be enough to conclude.

more

The following is a way to get the limit, dialog with the sage interpreter:

sage: M = matrix( QQ, 3, 3, [[0,5,10],[1/10,0,0],[0,1/2,0]] )
sage: M.charpoly().factor()
(x - 1) * (x^2 + x + 1/2)
sage: F.<a> = QuadraticField( -1 )
sage: M = M.change_ring( F )
sage: M.charpoly().factor()
(x - 1) * (x - 1/2*a + 1/2) * (x + 1/2*a + 1/2)
sage: J, S = M.jordan_form( transformation=True )
sage: J
[           1|           0|           0]
[------------+------------+------------]
[           0| 1/2*a - 1/2|           0]
[------------+------------+------------]
[           0|           0|-1/2*a - 1/2]
sage: S
[             1              1              1]
[          1/10 -1/10*a - 1/10  1/10*a - 1/10]
[          1/20         1/10*a        -1/10*a]
sage: S * J * S.inverse() == M
True


Now it is clear that $J^n$ converges to the diagonal matrix with entries $1,0,0$. We have $SJ^nS^{-1}=M^n$. So the limit is:

sage: S * diagonal_matrix( F, [1,0,0] ) * S.inverse()
[ 2/5    4    4]
[1/25  2/5  2/5]
[1/50  1/5  1/5]


We can test numerically:

sage: (M^100).n()
[ 0.399999999999999   4.00000000000000   4.00000000000000]
[0.0400000000000000  0.399999999999999  0.400000000000000]
[0.0200000000000000  0.200000000000000  0.199999999999999]


Note: Of course, there is no need to work over the gaussian field, where the characteristic polynomial of $M$ splits, but the prints are better.

more