8
$\begingroup$

I wish to compute the Taylor series expansion of the following iteration method $x_{k+1}=x_k-f'(x_k)^{-1}f(x_k)$ up to four terms of error ($e_k=x_k-\alpha$). When this is a scalar iteration, I very simply write the following

ClearAll["Global`*"]
f[e_] := df*(e^1 + c2 e^2 + c3 e^3 + c4 e^4 + c5 e^5);
fe = f[e];
f1e = f'[e];
Series[f1e^-1, {e, 0, 4}]*df // Simplify
u = e - Series[fe/f1e, {e, 0, 4}] // Simplify

and obtain correct results $$f(x_k)^{-1}=[1-2 \text{c2} e_k+ \left(4 \text{c2}^2-3 \text{c3}\right)e_k^2$$ $$-4 (2 c2^3 - 3 c2 c3 + c4) e_k^3 + (16 c2^4 - 36 c2^2 c3 + 9 c3^2 + 16 c2 c4 - 5 c5) e_k^4]f'({\alpha})^{-1}.$$ But, how to do this for the multi-dimensional case. That is, when the coefficients $c2,c3,c4$, and even $df=f'({\alpha})^{-1}$ are all matrices (note that e.g. $c2=\frac{1}{2!}f'(\alpha)^{-1}f^{(2)}(\alpha)$). By hand, I obtain the following correct results: $$ f(x_k)^{-1}=\left(I-2c2e_k+(4c2^2-3c3)e_k^2+(6c3c2+6c2c3-8c2^3-4c4)e_k^3\\ + (8c4c2+9c3^2+8c2c4-5c5-12c3c2^2-12c2c3c2-12c2^2c3+16c2^4)e_k^4\right)f'(\alpha)^{-1} $$ and \begin{equation} e_{k+1}=x_{k+1}-\alpha=-c2e_k^2+(2c2^2-2c3)e_k^3+(-4c2^3+4c2c3+3c3c2-3c4)e_k^4+\mathcal{O}(e_k^5). \end{equation}

So, my question is how to handle this non-commutative multiplication inside Series[] (e.g. $c2c3$ is not equal to $c3c2$ in this case)? I also tried to apply $**$, but I failed. Any suggestions or tricks to solve this problem will be appreciated fully.

I also saw the following posts in MathematicaStackExchange, but they were not useful for this problem, How to make noncommutative multiplication agree with commutative multiplication and Noncommutative multiply- expand expression.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
7
$\begingroup$

One idea is to use the function MatrixD from my answer to question det simplification that takes scalar derivatives of functions of matrices:

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]
    ]
]

Here are your expressions:

f[e_] := df*(IdentityMatrix[d] e^1 + c2 e^2 + c3 e^3 + c4 e^4 + c5 e^5);
fe = f[e];
f1e = f'[e];

I added the IdentityMatrix and assumed that df is just a scalar. Then:

MatrixD[Inverse[f1e], e] /. e->0

-Inverse[df IdentityMatrix[d]].(2 c2 df).Inverse[df IdentityMatrix[d]]

In order to make progress, we need to include assumptions about the coefficients:

$Assumptions = (c2 | c3 | c4 | c5) ∈ Matrices[{d, d}] && df ∈ Complexes;

The first few derivatives of Inverse[f1e] tensor expanded using the above assumptions:

TensorExpand[MatrixD[df Inverse[f1e], e] /. e->0]
TensorExpand[MatrixD[df Inverse[f1e], e, e] /. e->0]
TensorExpand[MatrixD[df Inverse[f1e], e, e, e] /. e->0]

-2 c2

-6 c3 + 8 MatrixPower[c2, 2]

-24 c4 + 36 c2.c3 + 36 c3.c2 - 48 MatrixPower[c2, 3]

So, the series expansion is:

TensorExpand[
    (df Inverse[f1e] /. e->0) +
    (MatrixD[df Inverse[f1e], e] /. e->0) e +
    (MatrixD[df Inverse[f1e], e, e] /. e->0) e^2/2! + 
    (MatrixD[df Inverse[f1e], e, e, e] /. e->0) e^3/3! +
    (MatrixD[df Inverse[f1e], e, e, e, e] /. e->0) e^4/4!
] + O[e]^5 //TeXForm

$\operatorname{IdentityMatrix}[d]-2 \operatorname{c2} e+e^2 (4 \operatorname{MatrixPower}[\operatorname{c2},2]-3 \operatorname{c3})+e^3 (-8 \operatorname{MatrixPower}[\operatorname{c2},3]+6 \operatorname{c2}.\operatorname{c3}+6 \operatorname{c3}.\operatorname{c2}-4 \operatorname{c4})+e^4 (-12 \operatorname{c3}.\operatorname{MatrixPower}[\operatorname{c2},2]-12 \operatorname{MatrixPower}[\operatorname{c2},2].\operatorname{c3}+16 \operatorname{MatrixPower}[\operatorname{c2},4]+9 \operatorname{MatrixPower}[\operatorname{c3},2]-12 \operatorname{c2}.\operatorname{c3}.\operatorname{c2}+8 \operatorname{c2}.\operatorname{c4}+8 \operatorname{c4}.\operatorname{c2}-5 \operatorname{c5})+O\left(e^5\right)$

which is in basic agreement with your inverse computation. A similar computation can be done for Inverse[f1e] . fe. However, I will instead present an approach teaching Series to use MatrixD for matrix functions. The internal function used by Series is System`Private`InternalSeries, so we teach this internal function about Dot and Inverse:

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
]

Now, we can use Series directly:

Series[df Inverse[f1e], {e, 0, 4}] //TeXForm

$\operatorname{IdentityMatrix}[d]-2 \operatorname{c2} e+e^2 (4 \operatorname{MatrixPower}[\operatorname{c2},2]-3 \operatorname{c3})+e^3 (-8 \operatorname{MatrixPower}[\operatorname{c2},3]+6 \operatorname{c2}.\operatorname{c3}+6 \operatorname{c3}.\operatorname{c2}-4 \operatorname{c4})+e^4 (-12 \operatorname{c3}.\operatorname{MatrixPower}[\operatorname{c2},2]-12 \operatorname{MatrixPower}[\operatorname{c2},2].\operatorname{c3}+16 \operatorname{MatrixPower}[\operatorname{c2},4]+9 \operatorname{MatrixPower}[\operatorname{c3},2]-12 \operatorname{c2}.\operatorname{c3}.\operatorname{c2}+8 \operatorname{c2}.\operatorname{c4}+8 \operatorname{c4}.\operatorname{c2}-5 \operatorname{c5})+O\left(e^5\right)$

in agreement with the previous result. Finally:

Series[Inverse[f1e] . fe, {e, 0, 4}] - e IdentityMatrix[d] //TeXForm

$-\operatorname{c2} e^2-2 e^3 (\operatorname{c3}-\operatorname{MatrixPower}[\operatorname{c2},2])+e^4 (-4 \operatorname{MatrixPower}[\operatorname{c2},3]+4 \operatorname{c2}.\operatorname{c3}+3 \operatorname{c3}.\operatorname{c2}-3 \operatorname{c4})+O\left(e^5\right)$

in agreement with your result.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Wow, a very useful response. The step-by-step way to get the multi-dimensional Taylor expansion is fruitful. $\endgroup$
    – Faz
    Commented Nov 16, 2018 at 17:48
2
$\begingroup$

This can be calculated by using NCAlgebra and this answer another question. The function

NCSeries[f_, {x_, x0_, n_}] := Block[{h}, 
SetNonCommutative[h]; 
Plus @@ (Table[1/i!, {i, 0, n}]*NestList[NCDirectionalD[#, {x, h}] &, f, n]) /. x -> x0 /. h -> x]

will calculate a formal power series of rational noncommutative functions. Your example is mixed, in the sense that e seems to commute but the coefficients do not. I propose that you lift e into a noncommutative variable ee

SNC[c1, c2, c3, c4, c5, ee] 
f[e_] := df*(e^1 + c2 ** e^2 + c3 ** e^3 + c4 ** e^4 + c5 ** e^5);
fe = f[ee]
f1e = f'[ee]   (* this works as expected in this simple case only! *)

use NCSeries to calculate the formal power series

series = NCSeries[ee - inv[f1e] ** fe, {ee, 0, 4}]

which you can then project back by replacing ee by the commutative e:

SetCommutative[e];
Collect[NCExpand[series /. ee -> e], e]

which will give the result you're looking for.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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