How can I expand expressions like $$(1-g(x)D)\frac{1}{1-f(x)D},$$ where $f(x),g(x)$ are some functions and $D$ is the shift operator defined by $Dh(x)=h(x+a)D$.

What I'm trying to do is to formally expand $$\frac{1}{1-f(x)D}=1+f(x)D+(f(x)D)(f(x)D)+...$$ and than multiply this with $(1-g(x)D)$, collecting all $D$ operators on the right most.

Is this possible with Mathematica?

  • 1
    $\begingroup$ You can add latex code with use of $, but Mathematica code is prefered. $\endgroup$
    – Feyre
    Oct 27, 2016 at 9:02
  • $\begingroup$ Your definition seems to contain an error. You mean $Dh(x)=h(x+a)$, right? $\endgroup$
    – Jens
    Oct 27, 2016 at 15:42
  • $\begingroup$ @Jens No, the definition is correct. You can see it considering the expression D (f(x) g(x)). $\endgroup$ Oct 27, 2016 at 15:46
  • $\begingroup$ @Dr.WolfgangHintze Sorry, but I'll wait for the OP to clarify because I don't see this as obvious. $\endgroup$
    – Jens
    Oct 27, 2016 at 15:51
  • 1
    $\begingroup$ @glS Again: this should be clarified in the question. It would affect how I'd implement this as code, and I don't want to guess. $\endgroup$
    – Jens
    Oct 27, 2016 at 17:36

4 Answers 4



This solution is a simpler altervative to that of gIS. It makes use only of the standard functions of Mathematica non-commutative multiplication (NCM) and replacement. Due to the uncommon features of the NCM, some care must be taken with linear combinations

In the first part we study the expression $\frac{1}{1-D f}$, which requires the powers of $(f D)$.

The second part is then devoted to the more general expression of the OP.

Part 1

We use the standard operation NonCommutativeMultiply[] as d and f do not commute (we write d instead of D to comply with Mathematica rules).

The shift operator will be implemented as the following replacement

r = d ** f[u_] -> f[u + a] ** d;

Now we have for the first few powers (notice that we have to use ReplaceRepeated[])

(f[x] ** d) //. r

(* Out[1111]= f[x] ** d *)

(f[x] ** d) ** (f[x] ** d) //. r

(* Out[1112]= f[x] ** f[a + x] ** d ** d *)

(f[x] ** d) ** (f[x] ** d) ** (f[x] ** d) //. r

(* Out[1113]= f[x] ** f[a + x] ** f[2 a + x] ** d ** d ** d *)

So the shift operator does what it should do. We can call the product where all d's are pushed through to the right "normal". Hence we know how to generate the normal product of the powers of (f[x] d)

Now the general power can be generated as

p[n_] := NonCommutativeMultiply @@ Table[f[x] ** d, {n}]

and the normal product is given by

pn[n_] := p[n] //. r


(* Out[1114]= f[x] ** f[a + x] ** f[2 a + x] ** d ** d ** d *)

We can generate the normal product of any power of the form (f[x] d). This completes the first part.

Part 2

The second part is not difficult. The only additional expression is the product g[x] d (f d)^n

For example (n=3)

q = g[x] ** d ** pn[3] //. r

(* Out[1138]= g[x] ** f[a + x] ** f[2 a + x] ** f[3 a + x] ** d ** d ** d ** d *)

To finalize all expressions we use this second replacement

rf = {d -> 1, NonCommutativeMultiply -> Times};

For example

q //. rf

(* Out[1139]= f[a + x] f[2 a + x] f[3 a + x] g[x] *)


q1 = pn[2] + pn[3]

(* Out[1140]= f[x] ** f[a + x] ** d ** d + f[x] ** f[a + x] ** f[2 a + x] ** d ** d ** d *)

q1 /. rf

(* Out[1141]= f[x] f[a + x] + f[x] f[a + x] f[2 a + x] *)

Some care has still to be taken in linear combinations: we need to apply the function Distribute[] and have to take (-1) as an expression to appear as a factor of the NCM.

The complete expression including g is then (in "finalized" form)

gf[n_] := Distribute[(1 + (-1) ** g[x] ** d) ** pn[n]] //. r //. rf



(* Out[1187]= 
f[x] f[a + x] f[2 a + x] - f[a + x] f[2 a + x] f[3 a + x] g[x] *)


1) The expression

$$\text{ff}=\frac{1}{1-D f}$$

can be written explicitly as

ff := 1 + Sum[Product[f[x + k a], {k, 0, n - 1}], {n, 1, \[Infinity]}]

$$\text{ff}\text{=}\sum _{n=1}^{\infty } \prod _{k=0}^{n-1} f(a k+x)+1$$

This expression can then be studied for further simplification depending on the function f.

2) In a comment, Jens pointed out that the term "shift operator" is reserved in standard literature as e.g. in quantum mechanics text books, and defined there as D f(x) = f(x+a) instead of D f(x) = f(x+a) D as in the OP.

Considering a typical expression of the OP, w = (D f)(D f), we see the difference

Standard use:

(D f)(D f) = f(x+a) D f = f(x+a) f(x+a) = f(x+a)^2  


(D f)(D f) = f(x+a) D D f = f(x+a) D f(x+a) D = f(x+a) f(x+2a) D^2

I have adopted the understanding of the OP. In standard use the problem is trivial.

3) Example




the exponential operator gives

$$\text{fe}=\exp (f(x) d)=(1-a)^{-\frac{x}{a}}$$


To implement something like this you first of all need to use some kind of noncommutative product, like NonCommutativeMultiply, otherwise terms get reordered automatically. Given that you also want to implement some special rules you can define a wrapper (which I called NCP as in Non Commutative Product) in which the symbols are handled in a special way. When the reordering is done you can replace NCP with a regular NonCommutativeMultiply, or just leave with NCP, depending on what you want to do with the final expression.

In the following code I set up things in such a way that only the symbols listed in targetFunctions are targeted by the shift operator, which I denoted wiht \[ScriptCapitalD]

(* targetFunction defined what are the symbols on which the shift operator should operate *)
ClearAll[targetFunctions, f, g, h];
targetFunctions = (f | g | h);
(* scalar is a pattern defining what we should consider as a "scalar" when handling the noncommutative product (for now it's just numeric values but you can add whatever you want) *)
scalar = (_?NumericQ);
(* powers of noncommutative products are recasted as flat noncommutative products *)
NCP /: Power[NCP[x__], n_Integer] := Apply[NCP][Join @@ ConstantArray[{x}, n]];
(* implement the behaviour we expect from a noncommutative product *)
NCP[] = 1;
NCP[left___, x : scalar, right___] := x NCP[left, right];
NCP[ll___, Times[x : scalar, f_], rr___] := x NCP[ll, f, rr];
NCP[left___, NCP[x__], right___] := NCP[left, x, right];
NCP[x_] := x;

ClearAll[applyShift, \[ScriptCapitalD]];
(* distribute NCP over sum *)
applyShift[HoldPattern[Plus[expr__]]] := Apply[Plus][applyShift /@ {expr}];
(* take out scalars *)
applyShift[Times[x : scalar, f_]] := x applyShift[f];
(* rule defining the action of the shift operator on the target functions *)
applyShift[NCP[left___, \[ScriptCapitalD][a_], (f : targetFunctions)[x_], 
    right___]] := applyShift @ NCP[left, f[x + a], \[ScriptCapitalD][a], right];
(* remove the special wrapper (delete this line if you want to keep NCP *)
applyShift[x___] := x /. NCP -> NonCommutativeMultiply;

and to see it in action:

(* formally expand geometric series and replace x with f[x] D[a] ... *)
Normal[Series[1/(1 - x), {x, 0, 4}]] /. {x -> NCP[f[x], \[ScriptCapitalD][a]]}
(* ... and apply shift rules *)
Output[1] = 1 + NCP[f[x], \[ScriptCapitalD][a]] + 
 NCP[f[x], \[ScriptCapitalD][a], f[x], \[ScriptCapitalD][a]] + 
 NCP[f[x], \[ScriptCapitalD][a], f[x], \[ScriptCapitalD][a], 
  f[x], \[ScriptCapitalD][a]] + 
 NCP[f[x], \[ScriptCapitalD][a], f[x], \[ScriptCapitalD][a], 
  f[x], \[ScriptCapitalD][a], f[x], \[ScriptCapitalD][a]]

Output[2] = 1 + f[x] ** \[ScriptCapitalD][a] + 
 f[x] ** f[a + x] ** \[ScriptCapitalD][a] ** \[ScriptCapitalD][a] + 
 f[x] ** f[a + x] ** 
  f[2 a + x] ** \[ScriptCapitalD][a] ** \[ScriptCapitalD][
   a] ** \[ScriptCapitalD][a] + 
 f[x] ** f[a + x] ** f[2 a + x] ** 
  f[3 a + x] ** \[ScriptCapitalD][a] ** \[ScriptCapitalD][
   a] ** \[ScriptCapitalD][a] ** \[ScriptCapitalD][a]

enter image description here


Taking on the idea of gIS and Dr. Wolfgang Hintze to lift the shift operator to a noncommutative operator, this can be done by using NCAlgebra and this answer another question. The function

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

will calculate the first n terms of the formal power series of a rational noncommutative function f around x0. Along with the rule

rule = {d ** f_[u_] -> f[u + a] ** d}

one can calculate terms of the noncommutative power series expansion, for example for degree 3:

exp = NCReplaceRepeated[NCSeries[(1 - g[x] ** d) ** inv[1 - f[x] ** d], {d, 0, 3}], rule];

which can get projected back into the commutative world as

BeginCommuteEverything[{d, f, g, a, x}];
Collect[exp, d]
(* 1 + d (f[x] - g[x]) + d^2 (f[x] f[a + x] - f[a + x] g[x]) + d^3 (f[x] f[a + x] f[2 a + x] - f[a + x] f[2 a + x] g[x]) *)

May be you should write expansion for fraction and then make a substitution

Series[1/(1-alpha), {alpha, 0, 5}]/.{alpha-> f[x]d}
  • $\begingroup$ But I want collect all $D$ to the right most at a give order, how could I find the coefficient? $\endgroup$
    – HH. Chen
    Oct 27, 2016 at 10:13
  • $\begingroup$ @HH.Chen If you want to collect all 'D' to the right then you should rewrite your expression. Now this expression is read as you first act with 'D', then multiply by 'f(x)' then again you act with 'D' and then multiply by 'g(x)'. It is better to rewrite your expression without two entries of 'D' $\endgroup$ Oct 27, 2016 at 14:46

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.