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!
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]
Really good solution! That's exactly what I was looking for!
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.
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.
Asked: 2017-12-20 09:17:31 +0100
Seen: 1,330 times
Last updated: Dec 28 '17
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
Is this homework ?
Yes, in the past it was! Now I'm using it as an exercise to learn a bit of Sage
Thank you so much for your answers!
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:
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).
Related: