# Replacing values of a function

I have recently starting using Mathematica and have recently come to what seems as an impasse in my understanding of the language. If this is too "tell me how to do it" I would certainly understand.

I have written a program that applies commutator logic to a sequence of creation and annihilation operators to produce a string of Kronecker deltas. I then need to apply the Kronecker delta restrictions to my operators. This is purely symbolic algebra, no real numbers need be produced at the end of the day.

For example,

(* KD[a,b] is a Kronecker delta with indices a and b. *)
(* KD Attributes: Orderless,NumericFunction *)
(* A[{x},{y}] is an operator with creation indices x and annihilation indices y*)

In:
A[{x},{y}] (KD[x,a1]*KD[y,a2]-KD[x,b1]*KD[y,b2])

Out:
A[{a1},{a2}]-A[{b1},{b2}]


What I have so far:

In:
rules=torules[(KD[x,a1]*KD[y,a2]]
A[{x},{y}]/.rules

Out:
{x->a1,y->a2}
A[{a1},{a2}]


This, however, is not quite what I need and converting from KD's to the rules list seems overly complicated. The main issues with this is that I cannot generate operator expressions that have plus within them. The structure of KD is not easily changed, but the input of the operator expressions is quite flexible.

Another take on this, I am not sure if it will help, but its my current direction.

In:
SetAttributes[A, NumericFunciton]
A[{x},{y}] (KD[x,a1]*KD[y,a2]-KD[x,b1]*KD[y,b2])

Out:
A[{x},{y}]*KD[x,a1]*KD[y,a2]- A[{x},{y}]*KD[x,b1]*KD[y,b2]


This should be a fairly simple replacement pattern, but I cannot get any further.

-
(1) If your KD is Orderless then how do you enforce that, say, A[{x},{a2}] KD[x,a1]*KD[y,a2] will evaluate to A[{a1},{a2}] rather than A[{x},{y}] or A[{x},{a2}] or A[{a1},{y}]? (2) What is your torules[] definition? –  Daniel Lichtblau Feb 6 '13 at 22:20
Have you thought about using one of the built-in deltas? Say, replace KD[x,a] with DiracDelta[x-a] and then just Integrate over {x,-Infinity,Infinity}. –  Xerxes Feb 6 '13 at 22:23
@DanielLichtblau (1)The deltas when applied to an operator in this way is similar to saying if x replace with a1 or if a1 replace with x in the example KD[x,a1]. The formalism is designed so that you will not have these issues. (2)It does not exist, I was trying to supply at least something that I thought of. –  Ophion Feb 6 '13 at 22:33
@Xerxes My understanding points towards these functions being used to compute numerical expressions within Mathematica. I am trying to use Mathematica to generate symbolic expressions that will be used by another program. –  Ophion Feb 6 '13 at 22:35
Integrate is symbolic. Perhaps you are thinking of NIntegrate? Try evaluating this to see what I mean: Integrate[ A[x, y] DiracDelta[x - a] DiracDelta[y - b], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, Assumptions -> {a, b} [Element] Reals] –  Xerxes Feb 6 '13 at 22:44

Here are several rules that can help to implement the Einstein convention:

deltaSimplifyRules =
{
(a_ + b_)*x_ :> (a*x + b*x),
KD[i_Symbol , j_]*
x_ /; (MemberQ[Attributes[i], Temporary] &&
MemberQ[x, i, Infinity]) :> (x /. i :> j),
KD[i_ , j_Symbol]*
x_ /; (MemberQ[Attributes[j], Temporary] &&
MemberQ[x, j, Infinity]) :> (x /. j :> i),
KD[i_Symbol , i_] :>  3 /; MemberQ[Attributes[i], Temporary],
KD[i_Symbol , j_]^2 /; MemberQ[Attributes[i], Temporary] :> KD[j, j],
KD[x__] :> Sort[KD[x]]
};

simplifyKD[x_] := x //. deltaSimplifyRules


now you can use simplifyKD to simplify your expression:

simplifyKD[
Module[{x, y},
A[{x}, {y}] (KD[x, a1]*KD[y, a2] - KD[x, b1]*KD[y, b2])  ]]

A[{a1}, {a2}] - A[{b1}, {b2}]


Module is used to create a temporary variable. Any temporary variable is summed over by these rules.

-
This works really well with some tweaking! One thing that I am confused about is if I set the KD string to a variable and then: simplifyKD[Module[{x, y}, A[{x}, {y}]*(KDstring)] This does not work. Its easy to work around, but interesting that its a problem. –  Ophion Feb 7 '13 at 14:23

Assuming you want to normal-order products of creation and annihilation operators, followed by tidying up the Kronecker deltas and the ordering of indices, you could do something like this (using an obvious notation):

Clear[opprod, del];


Define the basic commutation rule for moving annihilation operators to the right.

opprod[u___, an[i_], cr[j_], v___] := opprod[u, cr[j], an[i], v] + del[i, j] opprod[u, v];


Tidy up a corner case.

opprod[] = id;


Define rules for simplifying and tidying Kronecker deltas - there are other possibilities that could be used here.

del[i_, i_] = 1;
del[i_, j_] /; ! OrderedQ[{i, j}] := del[j, i];
del /: del[i_, j_]^_ := del[i, j];
del /: del[i_, j_] del[j_, k_] := del[i, j] del[i, k];
del /: del[i_, k_] del[j_, k_] := del[i, j] del[i, k];


Define a rule for standardising the indices between Kronecker deltas and operator products.

opprod /: del[i_, j_] opprod[u___, (x : (an | cr))[j_], v___] :=
del[i, j] opprod[u, x[i], v];


Define a rule for tidying up indices in operator products - this puts the indices of commuting operators in a canonical order.

opprod[u___, (x : (an | cr))[i_], (x : (an | cr))[j_], v___] /; !OrderedQ[{i, j}] :=
opprod[u, x[j], x[i], v];


Now try out these definitions on a non-trivial operator product.

opprod[an[k], an[j], cr[i], an[i], cr[j], cr[k]] // Expand

(* id del[i, j] + id del[i, k] + 2 id del[i, j] del[i, k] +
opprod[cr[i], an[i]] + 2 del[i, j] opprod[cr[i], an[i]] +
2 del[i, k] opprod[cr[i], an[i]] +
6 del[i, j] del[i, k] opprod[cr[i], an[i]] +
del[j, k] opprod[cr[i], an[i]] + del[i, k] opprod[cr[j], an[j]] +
del[i, j] opprod[cr[k], an[k]] +
opprod[cr[i], cr[j], an[i], an[j]] +
2 del[i, k] opprod[cr[i], cr[j], an[i], an[j]] +
2 del[j, k] opprod[cr[i], cr[j], an[i], an[j]] +
opprod[cr[i], cr[k], an[i], an[k]] +
2 del[i, j] opprod[cr[i], cr[k], an[i], an[k]] +
opprod[cr[i], cr[j], cr[k], an[i], an[j], an[k]] *)


I have quoted this result, not because it is necessarily interesting, but so you can verify that the above definitions are working as I had intended. I have cross-checked that this result is the same as what I get using another home-grown method that I have "up my sleeve". If it turns out to be incorrect then please tell me!

It is also very useful to use the Notation package to format the output in a human-readable way. For instance, you could try this (it looks much less messy when you enter it using the Notation package's palette):

<< Notation

Notation[ParsedBoxWrapper[
SubscriptBox["C",
"i_"]]\[DoubleLongLeftRightArrow]ParsedBoxWrapper[
RowBox[{"cr", "[", "i_", "]"}]]]
Notation[ParsedBoxWrapper[
SubscriptBox["A",
"i_"]]\[DoubleLongLeftRightArrow]ParsedBoxWrapper[
RowBox[{"an", "[", "i_", "]"}]]]
InfixNotation[ParsedBoxWrapper["\[SmallCircle]"], opprod]
Notation[ParsedBoxWrapper[
SubscriptBox["\[Delta]",
RowBox[{"i_", ",",
"j_"}]]]\[DoubleLongLeftRightArrow]ParsedBoxWrapper[
RowBox[{"del", "[", RowBox[{"i_", ",", "j_"}], "]"}]]]

-
This is pretty close to what I came up with. The /: and Notation is a neat trick and will definitely clean up the code some. I am curious though why del /: del[i_, j_]^_ := del[i, j]; is necessary? –  Ophion Feb 7 '13 at 14:31
When you use the commutator rule two or more times in succession it generates products of Kronecker deltas, so terms such as del[i,j]^2 (etc) can arise. The del /: del[i_, j_]^_ := del[i, j] definition simplifies these. –  Stephen Luttrell Feb 7 '13 at 17:17
Ah I was reading that as del /: del[i_, j_] ^= del[i, j]. One additional question, is there any benefit to the above over something like del[i_, j_]^_ := del[i, j] without the del /:? –  Ophion Feb 7 '13 at 17:37
The effect of del /: is to make the definition associate itself with del, which is what you would want. If you omit del /: then it defaults to trying to associate itself with the Head of del[i_, j_]^_ which is Power, which is not what you you would want. Try a_^b_ := f[a,b]` to see what happens. –  Stephen Luttrell Feb 8 '13 at 0:19