Formal expansion of matrices using Maple or Mathematica

Question

I'd like to evaluate some powers of sum of matrices, say I would like to evaluate $\left(A+B\right)^{n}$ with $A$ and $B$ some matrices. Since I'd like to go to high order, I'd like to use Maple or Mathematica to do that, but is there a way ?

To be more precise, how to calculate non-commutative product of commutative sum of objects in Maple or Mathematica ?

Of course, I want the formal expansion, for instance I want to calculate $\left(A+B\right)^{2}$ and get $$\left(A+B\right)\cdot\left(A+B\right)=A\cdot A+A\cdot B+B\cdot A+B\cdot B$$ at this level I do not need software, but at high order it would be convenient to have some. At no point I want to define what are $A$ or $B$.

Answer 1

In Maple you could do something like this below, though you might find that it's also possible using the Physics package.

You can of course adjust and extend this kind of thing. Perhaps it'll get you started.

Maple's expand command already knows to treat the infix &* operator as noncommuting. We can create a short procedure to turn positive powers into explicit products involving that &* operator.

restart:

ncexpand:=proc(ee)                                                  
  subsindets(ee, `^`(anything,posint),
             proc(u)
               local i;
               `&*`(seq(op(1,u),i=1..op(2,u)));
             end proc);
end proc:

ncexpand( (A+B)^2 );

                        (A + B) &* (A + B)

expand( ncexpand( (A+B)^2 ) );

             (A &* A) + (A &* B) + (B &* A) + (B &* B)

ncexpand( (A+B)^3 );

                    (A + B) &* (A + B) &* (A + B)

ncexpand( (A+B)^3 + 7*(A+B)^2 );

        ((A + B) &* (A + B) &* (A + B)) + 7 ((A + B) &* (A + B))

S:=expand(ncexpand( (A+B)^3 + 7*(A+B)^2 )):

Let's check that:

eval( eval(S, [A=Matrix([[1,2],[3,4]]),
               B=Matrix([[1,1],[2,2]])]),
      `&*`=`.` );

                            [291    369]
                            [          ]
                            [615    783]

eval( (A+B)^3 + 7*(A+B)^2,
      [A=Matrix([[1,2],[3,4]]),
       B=Matrix([[1,1],[2,2]])] );

                            [291    369]
                            [          ]
                            [615    783]