Ask Your Question
2

A^n, n=var(’n')

asked 2011-08-18 15:51:31 -0500

alreiff gravatar image

We have a matrix, for instance A = Matrix([[1,2],[3,4]]) How caculate A^n with n a variable (i.e symbolic calculus) and not n=50 or n = 100

edit retag flag offensive close merge delete

4 answers

Sort by » oldest newest most voted
1

answered 2012-06-25 04:17:26 -0500

pang gravatar image

If the matrix is diagonalizable, the following works:

sage: def pow(D,k):
...       return diagonal_matrix([l^k for l in D.diagonal()])

sage: D,P=M.jordan_form(transformation=True)
sage: show(P*pow(D,k)*~P)

It gets more messy if the matrix admits a jordan form in Q, and it needs some patching if some eigenvalues are in QQbar\QQ, because there are no canonical coercions between QQbar and SR, the likely parent of the final result. Code like the following did the trick:

sage: D,P = M.base_extend(QQbar).jordan_form(transformation = True)
sage: matriz = (P.apply_map(lambda l:SR(l))*powJordan(D,k).apply_map(lambda l:SR(l))*(~P).apply_map(lambda l:SR(l)))

apply_map is used because base_extend does not work, for the reasons mentioned above.

powJordan could be for example like this:

sage: def powJordan(D, k):
...       '''Eleva una matriz de Jordan (no necesariamente diagonal) a una variable simbolica'''
...       
...       dim = D.ncols()
...       M = matrix(SR, dim, dim)
...   
...       i = 0
...       while i<=dim-1:
...           #Autovalores simples, o con multiplicidad algebraica == multiplicidad geometrica
...           if (i== dim-1) or D[i,i+1] == 0: 
...               M[i,i] = SR(D[i,i])^k
...               i = i + 1
...           else: #buscamos la dimension del subespacio invariante
...               counter = 1
...               for j in range(i, dim-1):
...                   if D[j,j+1] == 0:
...                       break
...                   counter += 1
...               for j in srange(i, i+counter):
...                   for t in srange(j,(i+counter)):
...                       M[j,t] = binomial(k, t-j)*SR(D[i,i])^(k-(t-j))
...               i += counter
...           
...       return M

I have to say, I'm not fully satisfied with the solution, but it might work for you as it did for me.

edit flag offensive delete link more
1

answered 2011-08-18 16:10:48 -0500

kcrisman gravatar image

As you saw, this gives a NotImplementedError in Sage, so you can't do it directly.

This is possible in Maxima, though it just gives you a symbolic answer.

sage: maxima_console()
<snip>
(%i7) A: matrix([1,2],[3,4]);
                                   [ 1  2 ]
(%o7)                              [      ]
                                   [ 3  4 ]
(%i11) A.A;
                                  [ 7   10 ]
(%o11)                            [        ]
                                  [ 15  22 ]
(%i13) A^^2;
                                  [ 7   10 ]
(%o13)                            [        ]
                                  [ 15  22 ]
(%i14) A^^n;
                                          <n>
                                  [ 1  2 ]
(%o14)                            [      ]
                                  [ 3  4 ]

Unfortunately, I can't find this exponential command for matrices with letters, so we could use it by tab-completion in Sage and then use

sage: A = matrix([[1]])
sage: M = A._maxima_()
sage: M.command_name('n')

or something.

edit flag offensive delete link more
0

answered 2017-04-28 14:02:34 -0500

mforets gravatar image

From Sage v.8.0.beta3,

sage: A = Matrix([[1,2],[3,4]])
sage: n = SR.var('n')  
sage: A^n
[ 1/22*(-1/2*sqrt(33) + 5/2)^n*(sqrt(33) + 11) - 1/22*(1/2*sqrt(33) + 5/2)^n*(sqrt(33) - 11)                                                          2/33*sqrt(33)*(1/2*sqrt(33) + 5/2)^n - 2/33*sqrt(33)*(-1/2*sqrt(33) + 5/2)^n]
[-1/88*(-1/2*sqrt(33) + 5/2)^n*(sqrt(33) + 11)*(sqrt(33) - 3) - 1/88*(1/2*sqrt(33) + 5/2)^n*(sqrt(33) + 3)*(sqrt(33) - 11)                    1/66*sqrt(33)*(1/2*sqrt(33) + 5/2)^n*(sqrt(33) + 3) + 1/66*sqrt(33)*(-1/2*sqrt(33) + 5/2)^n*(sqrt(33) - 3)]

which is

$$ \left(\begin{array}{rr} \frac{1}{22} {\left(-\frac{1}{2} \sqrt{33} + \frac{5}{2}\right)}^{n} {\left(\sqrt{33} + 11\right)} - \frac{1}{22} {\left(\frac{1}{2} \sqrt{33} + \frac{5}{2}\right)}^{n} {\left(\sqrt{33} - 11\right)} & \frac{2}{33} \sqrt{33} {\left(\frac{1}{2} \sqrt{33} + \frac{5}{2}\right)}^{n} - \frac{2}{33} \sqrt{33} {\left(-\frac{1}{2} \sqrt{33} + \frac{5}{2}\right)}^{n} \\ -\frac{1}{88} {\left(-\frac{1}{2} \sqrt{33} + \frac{5}{2}\right)}^{n} {\left(\sqrt{33} + 11\right)} {\left(\sqrt{33} - 3\right)} - \frac{1}{88} {\left(\frac{1}{2} \sqrt{33} + \frac{5}{2}\right)}^{n} {\left(\sqrt{33} + 3\right)} {\left(\sqrt{33} - 11\right)} & \frac{1}{66} \sqrt{33} {\left(\frac{1}{2} \sqrt{33} + \frac{5}{2}\right)}^{n} {\left(\sqrt{33} + 3\right)} + \frac{1}{66} \sqrt{33} {\left(-\frac{1}{2} \sqrt{33} + \frac{5}{2}\right)}^{n} {\left(\sqrt{33} - 3\right)} \end{array}\right) $$

Also in the case of 2x2 matrices there is a derivation of the closed formula for $A^n$ in this note by K. S. Williams.

edit flag offensive delete link more
0

answered 2012-06-25 05:00:57 -0500

pang gravatar image

updated 2012-06-28 02:39:45 -0500

I've found another way that might work in some situations: use the exponential of the matrix:

$$e^M=\sum_{k\geq 0} \frac{1}{k!}M^k$$

then use derivatives to recover the powers of the matrix:

$$\left(\frac{\partial}{\partial k} \right)^k e^M = M^k$$

For example:

sage:  B = matrix(SR, 2, [1,2,3,4])
...    var('t')
...    f = exp(t*B)[1,1]
sage:  print n(f.derivative(t,4)(t=0))
634.000000000000
sage:  print n((B^4)[1,1])
634.000000000000

In order words, the k-th order derivative of the function f = exp(t*B)[1,1] at zero is the [1,1] coefficient of B^k. Maybe that solves your problem.

This could be useful in combination with LazyPowerSeries, as seen in the sage book in french.

Once you've got the function $f = (e^M)_{ij}$, you could get the function k -> [k-th order derivative of $f$] = [(i,j) entry of $M^k$] by Fourier transform, product with $\psi^k$, and inverse Fourier transform. But those integrals may be non-trivial...

edit flag offensive delete link more

Comments

Ah, pang! I will hijack this to ask what you think of the updates at #10637... they work for me...

kcrisman gravatar imagekcrisman ( 2012-06-25 05:56:38 -0500 )edit

Hi, karl: It's working now, I just need to doctest the full sage library and stuff like that. Currently my hardware is limited, so I cannot do this at any time.

pang gravatar imagepang ( 2012-06-28 02:34:55 -0500 )edit

I'll respond there.

kcrisman gravatar imagekcrisman ( 2012-06-28 04:08:59 -0500 )edit

Your Answer

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

Add Answer

Question Tools

Stats

Asked: 2011-08-18 15:51:31 -0500

Seen: 251 times

Last updated: Apr 28