Ask Your Question
1

Limit of a matrix power

asked 2017-12-20 02:17:31 -0600

dg.aragones gravatar image

updated 2017-12-20 17:08:44 -0600

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 flag offensive close merge delete

Comments

Is this homework ?

tmonteil gravatar imagetmonteil ( 2017-12-20 05:09:23 -0600 )edit

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

dg.aragones gravatar imagedg.aragones ( 2017-12-20 16:57:57 -0600 )edit

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:

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).

dg.aragones gravatar imagedg.aragones ( 2017-12-20 17:21:19 -0600 )edit

3 answers

Sort by » oldest newest most voted
3

answered 2017-12-26 16:06:32 -0600

slelievre gravatar image

updated 2017-12-27 20:13:53 -0600

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]
edit flag offensive delete link more

Comments

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

dg.aragones gravatar imagedg.aragones ( 2018-01-01 15:37:46 -0600 )edit
0

answered 2017-12-20 05:09:03 -0600

tmonteil gravatar image

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.

edit flag offensive delete link more
0

answered 2017-12-20 05:02:16 -0600

dan_fulea gravatar image

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.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-12-20 02:17:31 -0600

Seen: 75 times

Last updated: Dec 27 '17