How to implement the action of a shift operator on arbitrary functions?

Question

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?

Answer 1

Overview

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

pn[3]

(* 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] *)

or

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

Example

gf[3]

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

Discussion

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  

OP:

(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

Example

With

$$f(x)=x$$

the exponential operator gives

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

Answer 2

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) *)
ClearAll[scalar];
scalar = (_?NumericQ);
ClearAll[NCP];
(* 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 *)
applyShift@%
(*
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]
*)

Answer 3

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}, 
SetNonCommutative[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]) *)
EndCommuteEverything[]

Answer 4

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

Series[1/(1-alpha), {alpha, 0, 5}]/.{alpha-> f[x]d}