2
$\begingroup$

I am trying to obtain a power series expansion in some real parameter, but all the terms are arbitrary products of tensors. I.e. I want to expand an expression containing sums/products/powers of this parameter (say a) and arbitrary sums/products/powers of tensors as much as possible. For example

a((MatrixPower[X.Y + Y.X, 2] + a X.Y.Y.X))

should give

a(X.Y.X.Y + X.Y.Y.X + Y.X.X.Y + Y.X.Y.X) + a^2 X.Y.Y.X

After trying combinations of TensorExpand, TensorReduce, Distribute, Collect, Refine, Replace, .. I'm lost. For example I get a term like

$Assumptions = a \[Element] Reals && A1 \[Element] Matrices[{4, 4}] && A2 \[Element] Matrices[{4, 4}] && A3 \[Element] Matrices[{4, 4}] && A4 \[Element] Matrices[{4, 4}];

a^3 (A1.A2/2 + A1.A3/2 - A2.A1/2 + A2.A3/2 - A3.A1/2 - A3.A2/2
+ 1/12 A3.(a A1 + a A2 + 1/2 a^2 (A1.A2 - A2.A1) + 1/12 a^3 (A1.A1.A2
- 2 A1.A2.A1 + A2.A1.A1) + 1/12 a^3 (A1.A2.A2 - 2 A2.A1.A2 + A2.A2.A1)).
(a A1 + a A2 + 1/2 a^2 (A1.A2 - A2.A1) + 1/12 a^3 (A1.A1.A2 - 2 A1.A2.A1
+ A2.A1.A1) + 1/12 a^3 (A1.A2.A2 - 2 A2.A1.A2 + A2.A2.A1)))

after Collect, which is obviously not only a^3, but also a^4, a^5, a^6. I cannot get this into a form like

a^3(...) + a^4(...) + a^5(...) + a^6(...)

TensorExpand / TensorReduce do not do what I want because they write MatrixProduct[.., ..] whenever possible which still contains factors of a inside and then cannot be properly used for Collect. Distribute also works only in the simplest cases.

Is it somehow possible to do what I want?

$\endgroup$
2
$\begingroup$

You can use my code to implement series expansions of tensor objects in the question How to force Series[] to compute expansions by considering non commutative multiplication?, repeated at the end of this question. Your expression:

expr = a^3 (A1.A2/2 + A1.A3/2 - A2.A1/2 + A2.A3/2 - A3.A1/2 - A3.A2/2
+ 1/12 A3.(a A1 + a A2 + 1/2 a^2 (A1.A2 - A2.A1) + 1/12 a^3 (A1.A1.A2
- 2 A1.A2.A1 + A2.A1.A1) + 1/12 a^3 (A1.A2.A2 - 2 A2.A1.A2 + A2.A2.A1)).
(a A1 + a A2 + 1/2 a^2 (A1.A2 - A2.A1) + 1/12 a^3 (A1.A1.A2 - 2 A1.A2.A1
+ A2.A1.A1) + 1/12 a^3 (A1.A2.A2 - 2 A2.A1.A2 + A2.A2.A1)));

Using Series:

Series[expr, {a, 0, 6}] //TeXForm

$\frac{1}{2} a^3 (\operatorname{A1}.\operatorname{A2}-\operatorname{A2}.\operatorname{A1}+\operatorname{A1}. \operatorname{A3}-\operatorname{A3}.\operatorname{A1}+\operatorname{A2}.\operatorname{A3}-\operatorname{A3}.\operatorname{A2})+\frac{1}{12} a^5 (\operatorname{A3}.\operatorname{A1}.\operatorname{A2}+\operatorname{A3}.\operatorname{A2}. \operatorname{A1}+\operatorname{A3}.\operatorname{A1}.\operatorname{A1}+\operatorname{A3}.\operatorname{A2}.\operatorname{A2})+\frac{1}{24} a^6 (\operatorname{A3}.\operatorname{A1}.\operatorname{A1}.\operatorname{A2}+\operatorname{A3}. \operatorname{A1}.\operatorname{A2}.\operatorname{A2}-\operatorname{A3}.\operatorname{A2}.\operatorname{A1}.\operatorname{A1}-\operatorname{A3}.\operatorname{A2}.\operatorname{A2}.\operatorname{A1})+O\left(a^7\right)$

The code:

MatrixD[expr_, x__] := With[
    {old = OptionValue[SystemOptions[], "DifferentiationOptions"->"ExcludedFunctions"]},

    Internal`WithLocalSettings[
        SetSystemOptions["DifferentiationOptions"->"ExcludedFunctions"->Join[old, {Det, Inverse, Tr}]];

        Unprotect[D];
        (* handle list derivatives *)
        D[h:((Det|Tr|Inverse)[m_]), {z_, n_Integer}] := Nest[D[#, Replace[z, _List :> {z}]]&, h, n];
        D[h:((Det|Tr|Inverse)[m_]), {z_List}] := D[h, #]& /@ z;
        D[h:((Det|Tr|Inverse)[m_]), z_, y___] := D[D[h, z], y];

        (* define derivatives for Det, Tr, and Inverse *)
        D[Det[m_], z:Except[_List]] := Det[m] Tr[Inverse[m] . D[m,z]];
        D[Tr[m_], z:Except[_List]] := Tr[D[m,z]];
        D[Inverse[m_], z:Except[_List]] := -Inverse[m] . D[m, z] . Inverse[m],

        D[expr, x],

        SetSystemOptions["DifferentiationOptions"->"ExcludedFunctions"->old];
        Clear[D];
        Protect[D]
    ]
]

Unprotect[System`Private`InternalSeries];
System`Private`InternalSeries[a_Inverse|a_Dot, {e_,e0_,n_}] := SeriesData[
    e,
    e0,
    TensorReduce[TensorExpand[NestList[MatrixD[#,e]&,a,n]/.e->e0]]/Range[0,n]!,
    0,
    n+1,
    1
]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.