We need to understand how to simplify expressions in Yacas, and the best way is to try writing our own algebraic expression handler. In this chapter we shall consider a simple implementation of a particular non-commutative algebra called the Heisenberg algebra. This algebra was introduced by Dirac to develop quantum field theory. We won't explain any physics here, but instead we shall to delve somewhat deeper into the workings of Yacas.
We would like to be able to manipulate such expressions, expanding brackets, collecting similar terms and so on, while taking care to always keep the non-commuting terms in the correct order. For example, we want Yacas to automatically simplify 2**B(k1)**3**A(k2) to 6**B(k1)**A(k2). Our goal is not to implement a general package to tackle complicated non-commutative operations; we merely want to teach Yacas about these two kinds of "quantum objects" called A(k) and B(k), and we shall define one function that a physicist would need to apply to these objects. This function applied to any given expression containing A's and B's will compute something called a "vacuum expectation value", or "VEV" for short, of that expression. This function has "classical", i.e. commuting, values and is defined as follows: VEV of a commuting number is just that number, e.g. VEV(4)=4, VEV(delta(k-l))=delta(k-l); and VEV(X*A(k))=0, VEV(B(k)*X)=0 where X is any expression, commutative or not. It is straightforward to compute VEV of something that contains A's and B's: one just uses the commutation relations to move all B's to the left of all A's, and then applies the definition of VEV, simply throwing out any remaining q-numbers.
10 # A(_k)**B(_l) <-- B(l)**A(k) + delta(k-l); |
However, this is not going to work right away. In fact this will immediately give a syntax error because Yacas doesn't know yet about the new multiplication **. Let's fix that: we shall define a new infix operator with the same precedence as multiplication.
RuleBase("**", {x,y}); Infix("**", OpPrecedence("*")); |
Now we can use this new multiplication operator in expressions, and it doesn't evaluate to anything -- exactly what we need. But we find that things don't quite work:
In> A(_k)**B(_l) <-- B(l)**A(k)+delta(k-l); Out> True; In> A(x)**B(y) Out> B(l)**A(k)+delta(k-l); |
RuleBase("A", {k}); RuleBase("B", {k}); RuleBase("delta", {k}); |
In> A(y)**B(z)*2 Out> 2*(B(z)**A(y)+delta(y-z)); |
In> A(k)*2**B(l) Out> 2*A(k)**B(l); In> A(x)**A(y)**B(z) Out> A(x)**A(y)**B(z); In> (A(x)+B(x))**2**B(y)*3 Out> 3*(A(x)+B(x))**2**B(y); |
In> FullForm( A(k)*2**B(l) ) (** (* 2 (A k ))(B l )) Out> 2*A(k)**B(l); |
A solution for this problem is not to write rules for all possible cases (there are infinitely many cases) but to systematically reduce expressions to a canonical form. "Experience has shown that" (a phrase used when we can't come up with specific arguments) symbolic manipulation of unevaluated trees is not efficient unless these trees are forced to a pattern that reflects their semantics.
We should choose a canonical form for all such expressions in a way that makes our calculations -- namely, the function VEV() -- easier. In our case, our expressions contain two kinds of ingredients: normal, commutative numbers and maybe a number of noncommuting symbols A(k) and B(k) multiplied together with the "**" operator. It will not be possible to divide anything by A(k) or B(k) -- such division is undefined.
A possible canonical form for expressions with A's and B's is the following. All commutative numbers are moved to the left of the expression and grouped together as one factor; all non-commutative products are simplified to a single chain, all brackets expanded. A canonical expression should not contain any extra brackets in its non-commutative part. For example, (A(x)+B(x)*x)**B(y)*y**A(z) should be regrouped as a sum of two terms, (y)**(A(x)**(B(y))**A(z)) and (x*y)**(B(x)**(B(y))**A(z)). Here we wrote out all parentheses to show explicitly which operations are grouped. (We have chosen the grouping of non-commutative factors to go from left to right, however this does not seem to be an important choice.) On the screen this will look simply y ** A(x) ** B(y) and x*y** B(x) ** B(y) ** A(z) because we have defined the precedence of the "**" operator to be the same as that of the normal multiplication, so Yacas won't insert any more parentheses.
This canonical form will allow Yacas to apply all the usual rules on the commutative factor while cleanly separating all non-commutative parts for special treatment. Note that a commutative factor such as 2*x will be multiplied by a single non-commutative piece with "**".
The basic idea behind the "canonical form" is this: we should define our evaluation rules in such a way that any expression containing A(k) and B(k) will be always automatically reduced to the canonical form after one full evaluation. All functions on our new objects will assume that the object is already in the canonical form and should return objects in the same canonical form.
First, we need to distinguish "classical" terms from "quantum" ones. For this, we shall define a predicate IsQuantum() recursively, as follows:
/* Predicate IsQuantum(): will return True if the expression contains A(k) or B(k) and False otherwise */ 10 # IsQuantum(A(_x)) <-- True; 10 # IsQuantum(B(_x)) <-- True; /* Result of a binary operation may be Quantum */ 20 # IsQuantum(_x + _y) <-- IsQuantum(x) Or IsQuantum(y); 20 # IsQuantum(+ _y) <-- IsQuantum(y); 20 # IsQuantum(_x * _y) <-- IsQuantum(x) Or IsQuantum(y); 20 # IsQuantum(_x ** _y) <-- IsQuantum(x) Or IsQuantum(y); /* If none of the rules apply, the object is not Quantum */ 30 # IsQuantum(_x) <-- False; |
Now we shall construct rules that implement reduction to the canonical form. The rules will be given precedences, so that the reduction proceeds by clearly defined steps. All rules at a given precedence benefit from all simplifications at earlier precedences.
/* First, replace * by ** if one of the factors is Quantum to guard against user error */ 10 # (_x * _y)_(IsQuantum(x) Or IsQuantum(y)) <-- x ** y; /* Replace ** by * if neither of the factors is Quantum */ 10 # (_x ** _y)_(Not(IsQuantum(x) Or IsQuantum(y))) <-- x * y; /* Now we are guaranteed that ** is used between Quantum values */ /* Expand all brackets involving Quantum values */ 15 # (_x + _y) ** _z <-- x ** z + y ** z; 15 # _z ** (_x + _y) <-- z ** x + z ** y; /* Now we are guaranteed that there are no brackets next to "**" */ /* Regroup the ** multiplications toward the right */ 20 # (_x ** _y) ** _z <-- x ** (y ** z); /* Move classical factors to the left: first, inside brackets */ 30 # (x_IsQuantum ** _y)_(Not(IsQuantum(y))) <-- y ** x; /* Then, move across brackets: y and z are already ordered by the previous rule */ /* First, if we have Q ** (C ** Q) */ 35 # (x_IsQuantum ** (_y ** _z)) _(Not(IsQuantum(y))) <-- y ** (x ** z); /* Second, if we have C ** (C ** Q) */ 35 # (_x ** (_y ** _z))_(Not(IsQuantum(x) Or IsQuantum(y))) <-- (x*y) ** z; |
After we execute this in Yacas, all expressions involving additions and multiplications are automatically reduced to the canonical form. Extending these rules to subtractions and divisions is straightforward.
/* Commutation relation */ 40 # OrderBA(A(_k) ** B(_l)) <-- B(l)**A(k) + delta(k-l); 40 # OrderBA(A(_k) ** (B(_l) ** _x)) <-- OrderBA(OrderBA(A(k)**B(l)) ** x); /* Ordering simple terms */ 40 # OrderBA(_x)_(Not(IsQuantum(x))) <-- x; 40 # OrderBA(A(_k)) <-- A(k); 40 # OrderBA(B(_k)) <-- B(k); /* Sums of terms */ 40 # OrderBA(_x + _y) <-- OrderBA(x) + OrderBA(y); /* Product of a classical and a quantum value */ 40 # OrderBA(_x ** _y)_(Not(IsQuantum(x))) <-- x ** OrderBA(y); /* B() ** X : B is already at left, no need to order it */ 50 # OrderBA(B(_k) ** _x)<-- B(k) ** OrderBA(x); /* A() ** X : need to order X first */ 50 # OrderBA(A(_k) ** _x) <-- OrderBA(A(k) ** OrderBA(x)); |
These rules seem to be enough for our purposes. Note that the commutation relation is implemented by the first two rules; the first one is used by the second one which applies when interchanging factors A and B separated by brackets. This inconvenience of having to define several rules for what seems to be "one thing to do" is a consequence of tree-like structure of expressions in Yacas. It is perhaps the price we have to pay for conceptual simplicity of the design.
In> OrderBA(A(k)**A(l)) Error on line 1 in file [CommandLine] Line error occurred on: >>> Max evaluation stack depth reached. Please use MaxEvalDepth to increase the stack size as needed. |
50 # OrderBA(A(_k) ** _x) <-- OrderBAlate( A(k) ** OrderBA(x)); 55 # OrderBAlate(_x + _y) <-- OrderBAlate( x) + OrderBAlate(y); 55 # OrderBAlate(A(_k) ** B(_l)) <-- OrderBA(A(k)**B(l)); 55 # OrderBAlate(A(_k) ** (B(_l) ** _x)) <-- OrderBA(A(k)**(B(l)**x)); 60 # OrderBAlate(A(_k) ** _x) <-- A(k)**x; 65 # OrderBAlate(_x) <-- OrderBA(x); |
100 # VEV(_x) <-- VEVOrd(OrderBA(x)); /* Everything is expanded now, deal term by term */ 100 # VEVOrd(_x + _y) <-- VEVOrd(x) + VEVOrd(y); /* Now cancel all quantum terms */ 110 # VEVOrd(x_IsQuantum) <-- 0; /* Classical terms are left */ 120 # VEVOrd(_x) <-- x; |
Finally, we try some example calculations to test our rules:
In> OrderBA(A(x)*B(y)) Out> B(y)**A(x)+delta(x-y); In> OrderBA(A(x)*B(y)*B(z)) Out> B(y)**B(z)**A(x)+delta(x-z)**B(y) +delta(x-y)**B(z); In> VEV(A(k)*B(l)) Out> delta(k-l); In> VEV(A(k)*B(l)*A(x)*B(y)) Out> delta(k-l)*delta(x-y); In> VEV(A(k)*A(l)*B(x)*B(y)) Out> delta(l-y)*delta(k-x)+delta(l-x) *delta(k-y); |