Matrix multiplication in Maple

Question

I have defined my matrix M in Maple as follows:

M := Matrix(2, 2, [[-a, a], [b, -b]]);

Then I want to derive M^k in terms of M. The answer I am supposed to get is M^k = M*[-(a+b)]^{k-1}. When I enter the command:

M^k;

I am getting the following undesired output

Matrix(2, 2, [[-a, a], [b, -b]])^k

How to solve my problem? Picture from my Maple showing what I have tried.

Answer 1

Here is one way to get such a result:

M := Matrix(2, 2, [[-a, a], [b, -b]]):

simplify(LinearAlgebra:-MatrixPower(M,k));

      [           (k - 1)            (k - 1) ]
      [-a (-b - a)         a (-b - a)        ]
      [                                      ]
      [          (k - 1)              (k - 1)]
      [b (-b - a)          -b (-b - a)       ]

You can easily show that this equals M*(-b-a)^(k-1).

[edit] It occurs to me that this might be a coursework question, and as such you might be interested in a more step-by-step approach instead of just what result the MatrixPower command provides.

Let's diagonalize the Matrix M, checking results as we go.

restart;
with(LinearAlgebra):

M := Matrix([[-a, a], [b, -b]]);

               [-a  a ]
          M := [      ]
               [b   -b]

(evals,P) := Eigenvectors(M);

                          [  a   ]
                [-b - a]  [- -  1]
    evals, P := [      ], [  b   ]
                [  0   ]  [      ]
                          [ 1   1]

Diag := DiagonalMatrix(evals);

                [-b - a  0]
        Diag := [         ]
                [  0     0]

Check that P.Diag.P^(-1) equals M.

simplify( P . Diag . P^(-1) - M );

             [0  0]
             [    ]
             [0  0]

It turns out that using Adjoint(P)/Determinant(P) is slightly easier to use here than P^(-1). But first we'll check they're equivalent, which might be new to you.

Adjoint(P);

          [1   -1 ]
          [       ]
          [      a]
          [-1  - -]
          [      b]

detP := Determinant(P);

                 b + a
       detP := - -----
                   b  

simplify( Adjoint(P)/detP - P^(-1) );

             [0  0]
             [    ]
             [0  0]

We can also use denom(detP)/numer(detP) instead of 1/detP.

We want M^k. But that's the same as,

( P . Diag . P^(-1) )^k

Now, the "big idea:. If you multiply P.Diag.P^(-1) by itself k times then you get various instances of P^(-1).P in that product of terms. But P^(-1).P will just produce the identity matrix. So (P.Diag.P^(-1))^k will telescope down to,

P . Diag^k . P^(-1)

And here Diag^k is,

Diag^~k;

     [        k   ]
     [(-b - a)   0]
     [            ]
     [    0      0]

which is the same as,

Matrix([[1,0],[0,0]])*(-b-a)^k;

      [        k   ]
      [(-b - a)   0]
      [            ]
      [    0      0]

So M^k is the same as,

P . Matrix([[1,0],[0,0]])*(-b-a)^k . Adjoint(P)
  * denom(detP)/numer(detP);

   [            k             k ]
   [  a (-b - a)    a (-b - a)  ]
   [- -----------   ----------- ]
   [    -b - a        -b - a    ]
   [                            ]
   [           k               k]
   [ b (-b - a)      b (-b - a) ]
   [ -----------   - -----------]
   [   -b - a          -b - a   ]

You could check that that simplifies to the earlier MatrixPower result.

We could also rearrange the terms in that last product,

P . Matrix([[1,0],[0,0]]) . Adjoint(P) * denom(detP)
  * ((-b-a)^k) / numer(detP);

That last product can be made into two terms.

The "first term" turns out to be simply M,

P . Matrix([[1,0],[0,0]]) . Adjoint(P)
  * denom(detP);

         [-a  a ]
         [      ]
         [b   -b]

The "second term" is,

simplify( ((-b-a)^k) / numer(detP) );

              (k - 1)
      (-b - a)       

So M^k is the same as that first term times that second term, ie. M * (-b-a)^(k-1).