# Virasoro algebra in sympy

99 views

### glt

Oct 5, 2008, 7:54:52 AM10/5/08
to sympy
Hello all,

I am looking for a quick and short way to implement computations with
the Virasoro algebra in sympy. What i want to achieve is pretty basic:
i want to define symbolic noncommutative objects L(n), where n is an
integer, and have sympy handle products of L(n)'s as follows:

* On two-terms products L(n)L(m):
* if n<=m, do nothing
* otherwise, replace L(n)L(m) by L(m)L(n) + [L(n),L(m)] where [ ]
is the commutator of the Virasoro operators (expressed as another L
operator plus possibly another term)

* On products L(n)L(m)...L(p): recursively apply the two-terms rule,
so that n, m, .. p are in increasing order

So i started by subclassing Basic and ArithMeths in an L class, with
is_commutative=False. Now i'm stuck the the definition of the product
rule. I guess this should happen somewhere in Mul.flatten, but what is
the proper way of doing it? Subclassing Mul to define my own
multiplication? How can then i couple it to the L objects?

Another thing i have thought about would be representing the product
of L objects as a Basic object itself, with its own rules, but i'm not

How would you guys do it?

Tom

### glt

Oct 5, 2008, 7:54:52 AM10/5/08
to sympy

### glt

Oct 5, 2008, 7:54:52 AM10/5/08
to sympy

### Ondrej Certik

Oct 5, 2008, 9:31:21 AM10/5/08
Hi Tom!

On Sun, Oct 5, 2008 at 1:54 PM, glt <t.gu...@gmail.com> wrote:
>
> Hello all,
>
> I am looking for a quick and short way to implement computations with
> the Virasoro algebra in sympy. What i want to achieve is pretty basic:
> i want to define symbolic noncommutative objects L(n), where n is an
> integer, and have sympy handle products of L(n)'s as follows:
>
> * On two-terms products L(n)L(m):
> * if n<=m, do nothing
> * otherwise, replace L(n)L(m) by L(m)L(n) + [L(n),L(m)] where [ ]
> is the commutator of the Virasoro operators (expressed as another L
> operator plus possibly another term)
>
> * On products L(n)L(m)...L(p): recursively apply the two-terms rule,
> so that n, m, .. p are in increasing order

Thanks for doing it, that's very useful. What other operations would
you like to do with it?

Looking here, it's just a Lie algebra, right?

http://en.wikipedia.org/wiki/Virasoro_algebra

>
> So i started by subclassing Basic and ArithMeths in an L class, with
> is_commutative=False. Now i'm stuck the the definition of the product
> rule. I guess this should happen somewhere in Mul.flatten, but what is
> the proper way of doing it? Subclassing Mul to define my own
> multiplication? How can then i couple it to the L objects?
>
> Another thing i have thought about would be representing the product
> of L objects as a Basic object itself, with its own rules, but i'm not
>
> How would you guys do it?

I would suggest you either subclass Mul, or just "fork it", e.g.
create VirasoroMul by copying the current Mul and then customizing it
to do what you want. If you have problems, just create a function that
takes some set of L(n)L(m)...L(p) and returns the simplifications
using the algorithm you described. After you have some working
prototype, we'll help you integrate it well in sympy, with the patch
and getting it in.

In the past, we played with having a NCMul like ginac, then we
switched to just having one Mul and the commutative=False assumption,
currently attached to the symbols directly, in the future this will be
handled by the assumptions system, just like in Mathematica. If you
have a better idea how this could/should be handled, please share it.

Looking forward,
Ondrej

### Fredrik Johansson

Oct 5, 2008, 11:40:19 AM10/5/08
On Sun, Oct 5, 2008 at 3:31 PM, Ondrej Certik <ond...@certik.cz> wrote:
> In the past, we played with having a NCMul like ginac, then we
> switched to just having one Mul and the commutative=False assumption,
> currently attached to the symbols directly, in the future this will be
> handled by the assumptions system, just like in Mathematica. If you
> have a better idea how this could/should be handled, please share it.

But note that Mathematica handles noncommutative multiplication with a
separate "class", not with assumptions. I think a NCMul class would
be a better solution for SymPy too.

Fredrik

### glt

Oct 6, 2008, 3:39:40 PM10/6/08
to sympy
Hi, thank you all for your replies.

On Oct 5, 3:31 pm, "Ondrej Certik" <ond...@certik.cz> wrote:
>
> Thanks for doing it, that's very useful. What other operations would
> you like to do with it?

For now, essentially use sympy to sort the products properly,
expanding all the commutators, etc, which is a real pain to do by hand
for large products.

I would also be interested in defining anticommuting symbols, but I
guess that's really kind of the same problem.

>
> Looking here, it's just a Lie algebra, right?
>
> http://en.wikipedia.org/wiki/Virasoro_algebra
>
>

Yes it is. The idea of sorting the product terms by increasing order
is that it makes it easy to see which terms vanish on a lowest weight
representation vector of the algebra. You just have to discard all
products ending with an L_i with i>0.

>
> I would suggest you either subclass Mul, or just "fork it", e.g.
> create VirasoroMul by copying the current Mul and then customizing it
> to do what you want. If you have problems, just create a function that
> takes some set of L(n)L(m)...L(p) and returns the simplifications
> using the algorithm you described. After you have some working
> prototype, we'll help you integrate it well in sympy, with the patch
> and getting it in.

What i have done right now is very basic but works and does the job
(sort of). I've defined an L class which basically defines a product
L(a)L(b)...L(p)L(q), constructed in python as L(a,b,...,p,q).

Here are a few lines demonstrating how it works:

############

from sympy import *

class L(Basic):
is_commutative = False

def sort_virasoro(self, c):
"""Returns the product of Virasoro operators as a sum of
products of L operators in order of increasing indices"""
idx = self.args

# Find index of first pair of operators to swap
try:
iswap = [p>q for p,q in zip(idx[:-1],idx[1:])].index(True)
except:
return self

# Split L(a)L(b)...L(m)L(n)...L(p)L(q) into 4 parts:
# pre = L(a)L(b)..., L(m), L(n), and post = ...L(p)L(q)
pre = idx[:iswap]
post = idx[iswap+2:]
m = idx[iswap]
n = idx[iswap+1]

# This is the swapped term L(a)L(b)...L(n)L(m)...L(p)L(q)
out = L(*(pre + (n,m) + post)).sort_virasoro(c)

# This is the commutator L(a)L(b)...L(m+n)...L(p)L(q)*(m-n)
if m-n!=0:
out += (m-n)*L(*(pre+(m+n,)+post)).sort_virasoro(c)

# This is the central term
if m+n==0:
if len(pre+post)==0:
Lterm = 1
else:
Lterm = L(*(pre + post)).sort_virasoro(c)

out += Lterm * c/12*(m**3-m)

return out

if __name__=="__main__":
c = Symbol("c") # The central charge
foo = L(-7,3,-3,-4,-5).sort_virasoro(c)
print foo.expand()

############

Note that it's not optimal in many aspects (e.g. the sorting is
totally recursive and suboptimal) but it works.

I don't know how this kind of code could be integrated tighter with,
say Mul, so that you could just define L(3)*L(2)*L(1), and have sympy
automatically put it into this kind of "canonical" ordering for
example. Do you have insights on this? I've tried to look into the Mul
code but have not made it very far yet.

My education is much more in physics than in programming, but I'm
quite interested in sympy and its flexible use for this kind of
symbolic computation.

Tom

### Ondrej Certik

Oct 6, 2008, 4:02:25 PM10/6/08
On Mon, Oct 6, 2008 at 9:39 PM, glt <t.gu...@gmail.com> wrote:
>
> Hi, thank you all for your replies.
>
> On Oct 5, 3:31 pm, "Ondrej Certik" <ond...@certik.cz> wrote:
>>
>> Thanks for doing it, that's very useful. What other operations would
>> you like to do with it?
>
> For now, essentially use sympy to sort the products properly,
> expanding all the commutators, etc, which is a real pain to do by hand
> for large products.
>
> I would also be interested in defining anticommuting symbols, but I
> guess that's really kind of the same problem.

Right, thanks for the info.

>
>>
>> Looking here, it's just a Lie algebra, right?
>>
>> http://en.wikipedia.org/wiki/Virasoro_algebra
>>
>>
>
> Yes it is. The idea of sorting the product terms by increasing order
> is that it makes it easy to see which terms vanish on a lowest weight
> representation vector of the algebra. You just have to discard all
> products ending with an L_i with i>0.

I see. I had a similar problem with calculating cross-sections in QED
manually using creation and anihilation operators, one also needs to
employ commutation/anticommutation relations to simplify it. So SymPy
should be able to do this kind of manipulations.

>
>>
>> I would suggest you either subclass Mul, or just "fork it", e.g.
>> create VirasoroMul by copying the current Mul and then customizing it
>> to do what you want. If you have problems, just create a function that
>> takes some set of L(n)L(m)...L(p) and returns the simplifications
>> using the algorithm you described. After you have some working
>> prototype, we'll help you integrate it well in sympy, with the patch
>> and getting it in.
>
> What i have done right now is very basic but works and does the job
> (sort of). I've defined an L class which basically defines a product
> L(a)L(b)...L(p)L(q), constructed in python as L(a,b,...,p,q).
>
> Here are a few lines demonstrating how it works:
>
> ############
>
> from sympy import *
>
> class L(Basic):
> is_commutative = False
>
> def sort_virasoro(self, c):
> """Returns the product of Virasoro operators as a sum of
> products of L operators in order of increasing indices"""
> idx = self.args
>
> # Find index of first pair of operators to swap
> try:
> iswap = [p>q for p,q in zip(idx[:-1],idx[1:])].index(True)
> except:
> return self

^^^ Don't use empty "except:", because this catches all exceptions,
e.g. it hides problems. I guess you are catching the IndexError from
.index()?

>
> # Split L(a)L(b)...L(m)L(n)...L(p)L(q) into 4 parts:
> # pre = L(a)L(b)..., L(m), L(n), and post = ...L(p)L(q)
> pre = idx[:iswap]
> post = idx[iswap+2:]
> m = idx[iswap]
> n = idx[iswap+1]
>
> # This is the swapped term L(a)L(b)...L(n)L(m)...L(p)L(q)
> out = L(*(pre + (n,m) + post)).sort_virasoro(c)
>
> # This is the commutator L(a)L(b)...L(m+n)...L(p)L(q)*(m-n)
> if m-n!=0:
> out += (m-n)*L(*(pre+(m+n,)+post)).sort_virasoro(c)
>
> # This is the central term
> if m+n==0:
> if len(pre+post)==0:
> Lterm = 1
> else:
> Lterm = L(*(pre + post)).sort_virasoro(c)
>
> out += Lterm * c/12*(m**3-m)
>
> return out
>
> if __name__=="__main__":
> c = Symbol("c") # The central charge
> foo = L(-7,3,-3,-4,-5).sort_virasoro(c)
> print foo.expand()

Otherwise looks good.

>
> ############
>
>
> Note that it's not optimal in many aspects (e.g. the sorting is
> totally recursive and suboptimal) but it works.

Excellent, thanks for the code. I created an issue to get this into sympy:

>
> I don't know how this kind of code could be integrated tighter with,
> say Mul, so that you could just define L(3)*L(2)*L(1), and have sympy
> automatically put it into this kind of "canonical" ordering for
> example. Do you have insights on this? I've tried to look into the Mul
> code but have not made it very far yet.

Couple questions first:

1) should this happen automatically, or maybe only after calling some
function on the product?

If the latter, then all that is needed is to parse the noncummative
Mul expression, plug it into your function and construct the result
back. And integrate it nicely in sympy somewhere.

2) How is this implemented in other symbolic software, e.g.
Mathematica or Maple? (If you have experience with them.)

I learned it pays off to investigate other systems first, before
trying to come up with something ourselves.

> My education is much more in physics than in programming, but I'm
> quite interested in sympy and its flexible use for this kind of
> symbolic computation.

The same here, I also studied physics and not programming. What is
your application of this? E.g. what kind of stuff do you want to
calculate? My own little aim with SymPy is to be able to do basically
any calculation that one does in theoretical physics, but we still
have a long way to go. So I am glad you need this noncommutative
stuff, because it's just waiting to get polished and more widely used,
so if we can get SymPy do what you want, it will help us a lot.

Ondrej

### Brian Granger

Oct 10, 2008, 8:15:45 PM10/10/08
to sympy
Ondrej,

On a related note...

You mention that you have done QED related things with Sympy. Have
you ever implemented creation and annihilation operators along with
fock states in sympy?

Cheers,

Brian

On Oct 6, 1:02 pm, "Ondrej Certik" <ond...@certik.cz> wrote:

### Ondrej Certik

Oct 11, 2008, 5:49:08 AM10/11/08
On Sat, Oct 11, 2008 at 2:15 AM, Brian Granger <elliso...@gmail.com> wrote:
>
> Ondrej,
>
> On a related note...
>
> You mention that you have done QED related things with Sympy. Have
> you ever implemented creation and annihilation operators along with
> fock states in sympy?

Not yet, but I want to play with this in sympy. So far there is only:

\$ examples/qft.py

which calculates some simple scattering amplitude by direct
multiplications of the dirac matrices. Do you have some particular
example you'd like to play with? If so, I can try to implement that,
it should not be difficult.

Ondrej

### Brian Granger

Oct 11, 2008, 1:20:50 PM10/11/08
> Not yet, but I want to play with this in sympy. So far there is only:
>
> \$ examples/qft.py
>
> which calculates some simple scattering amplitude by direct
> multiplications of the dirac matrices. Do you have some particular
> example you'd like to play with? If so, I can try to implement that,
> it should not be difficult.

I would like to implement the bosonic and fermionic operators (and
associated many particle fock states) for the many-particle
(non-relativistic) Schrodinger equation. More specifically, like
those in Fetter and Walecka many body quantum text. As a start I
would be happy with the operators for a single state. This would
involve three basic entities:

a - annihilation operator
number state - |n>

along with either the commutation relations (bosons) or
anticommutation relations. I need to run right now, and I am just
starting to play with sympy, but there is a lot of potential for this
type of thing with sympy and it could have a significant impact.

Cheers,

Brian

> Ondrej
>
> >
>

### Ondrej Certik

Oct 12, 2008, 11:30:21 AM10/12/08