Consider the vector space \({\mathcal A}\) of linear operators on univariate functions; \({\mathcal A}\) can be made into an algebra where multiplication (denoted by juxtaposition) is defined as operator composition:

\[({\mathcal O}_1{\mathcal O}_2)f={\mathcal O}_1({\mathcal O}_2f)\]

We consider the algebra generated by by \(\partial\colon f\longrightarrow f'\) [that is, \((\partial f)(x) = f'(x)\)] and \(x\colon f\longrightarrow xf\) [that is, \((xf)(x) = xf(x)\)]. This is known as the (first) Weyl algebra. We observe that the Weyl algebra is not commutative: \(\partial xf=(xf)'=f+xf'\) but \(x\partial f=xf'\), so \(\partial x=x\partial+1\). The algebra generated by \(\left\lbrace x,\partial\right\rbrace\) will include elements such as \(7\partial + 4\partial x\partial^3 x\), which would map \(f\) to \(7f' + 4\left(x\left(xf\right)'''\right)'\).

It can be shown that any element of the Weyl algebra can be expressed in the standard form

\[ \sum_i a_i \partial^{p_i}x^{q_i}\]

for real \(a_i\) and nonnegative integers \(p_i,q_i\). Converting a general word to standard form is not straightforward but we have

\[ \partial x^n = x^n\partial + nx^{n-1}\]

and

\[ \partial^n x = x\partial^n + n\partial^{n-1}.\]

We can apply these rules recursively to find standard form for products \((\partial^i x^j)(\partial^l x^m)\). Alternatively we may follow Wolf and use the fact that

\[ (\partial^i x^j)(\partial^lx^m)= \sum_{r=0}^j{j\choose r}{l\choose r}\partial^{i+l-r}x^{j+m-r}.\]

All this is done in the package. Consider \(d_1=\partial x + 2\partial^3\) and \(d_2=3+7\partial -5x^2\partial^2\). Observing that \(d_1\) and \(d_2\) are in standard form, package idiom to create these operators would be:

`library(weyl)`

```
d1 <- weyl(spray(rbind(c(1,1),c(0,3)),c(1,2)))
d2 <- weyl(spray(rbind(c(0,0),c(1,0),c(2,2)),c(3,7,-5)))
d1
```

```
## A member of the Weyl algebra:
## x d val
## 0 3 = 2
## 1 1 = 1
```

`d2`

```
## A member of the Weyl algebra:
## x d val
## 2 2 = -5
## 1 0 = 7
## 0 0 = 3
```

Observe that, like the `spray`

package, the order of the matrix rows is not defined. We may apply the usual rules of arithmetic to these objects:

`d1*d2`

```
## A member of the Weyl algebra:
## x d val
## 2 5 = -10
## 1 4 = -60
## 2 2 = -10
## 0 3 = -54
## 1 3 = 14
## 1 0 = 7
## 0 2 = 42
## 3 3 = -5
## 2 1 = 7
## 1 1 = 3
```

We may choose to display results in symbolic form:

```
options(polyform=TRUE)
(d1^2 + d2) * (d2 - 3*d1)
```

```
## A member of the Weyl algebra:
## 9 -63*x^3*d^2 +168*d^5 +20*x^4*d^4 +87*x^3*d^3 -70*x^2*d -24*d^9
## +126*d^2 +49*x^2 -9*x*d -20*x^2*d^8 +112*x*d^3 -234*d^3 +49*x
## -178*x^2*d^5 -480*x*d^4 +24*x^2*d^2 -696*d^6 +28*x*d^6 -276*x*d^7
## -20*x^3*d^6 +28*x^2*d^4
```

Mathematica can deal with operators and we may compare the two systems’ results for \(\partial^2x\partial x^2\):

`In[1] := D[D[x*D[x^2*f[x],x],x],x] // Expand`

`Out[1] := 4 f[x] + 14 x f'[x] + 8 x^2 f''[x] + x^3f'''[x]`

```
x <- weyl(cbind(0,1))
D <- weyl(cbind(1,0))
x^2*D*x*D^2
```

```
## A member of the Weyl algebra:
## 4 +x^3*d^3 +8*x^2*d^2 +14*x*d
```

The package supports multivariate Weyl algebras. Consider:

```
options(polyform=FALSE) # revert to default print method
set.seed(0)
x <- rweyl()
x
```

```
## A member of the Weyl algebra:
## x y z dx dy dz val
## 2 0 1 2 0 1 = 3
## 0 1 2 2 0 1 = 2
## 1 0 2 1 0 1 = 1
```

Above, object `x`

is a member of the operator algebra generated by \(\left\lbrace\partial_x,\partial_y,\partial_z,x,y,z\right\rbrace\). We may verify associativity of multiplication:

```
x <- rweyl(n=1,d=2)
y <- rweyl(n=2,d=2)
z <- rweyl(n=3,d=2)
options(polyform=TRUE)
x*(y*z)
```

```
## A member of the Weyl algebra:
## +6*x*y^2*dx^2*dy +36*x*y^3*dx*dy^3 +2*x^2*y^4*dx*dy^4
## +4*x^2*y^2*dx*dy^5 +36*x*y^2*dx*dy^2 +x*y^4*dx*dy^3 +4*x*y^3*dx*dy^2
## +x^2*y^4*dx^2*dy^3 +12*x^2*y^2*dx*dy^2 +2*x^2*y^2*dx^2*dy^4
## +2*x*y^2*dx*dy^4 +2*x*y^2*dx*dy +2*x^2*y^2*dx^2*dy +4*x^2*y^3*dx^2*dy^2
## +3*x*y^4*dx^2*dy^3 +12*x^2*y^3*dx*dy^3 +6*x*y^4*dx*dy^4
## +12*x*y^3*dx^2*dy^2
```

`(x*y)*z`

```
## A member of the Weyl algebra:
## +2*x*y^2*dx*dy^4 +2*x^2*y^2*dx^2*dy^4 +3*x*y^4*dx^2*dy^3
## +x^2*y^4*dx^2*dy^3 +12*x^2*y^2*dx*dy^2 +x*y^4*dx*dy^3
## +2*x^2*y^4*dx*dy^4 +4*x^2*y^3*dx^2*dy^2 +2*x^2*y^2*dx^2*dy
## +2*x*y^2*dx*dy +4*x*y^3*dx*dy^2 +12*x^2*y^3*dx*dy^3 +6*x*y^4*dx*dy^4
## +12*x*y^3*dx^2*dy^2 +6*x*y^2*dx^2*dy +36*x*y^3*dx*dy^3
## +36*x*y^2*dx*dy^2 +4*x^2*y^2*dx*dy^5
```

Comparing the two results above, we see that they apparently differ. But the apparent difference is due to the fact that the terms appear in a different order, a feature that is not algebraically meaningful. We may verify that the expressions are indeed algebraically identical:

`x*(y*z) - (x*y)*z`

```
## A member of the Weyl algebra:
## the NULL multinomial of arity 4
```

A *derivation* \(D\) of an algebra \({\mathcal A}\) is a linear operator that satisfies \(D(d_1d_2)=d_1D(d_2) + D(d_1)d_2\), for every \(d_1,d_2\in{\mathcal A}\). If a derivation is of the form \(D(d) = [d,f] = df-fd\) for some fixed \(f\in{\mathcal A}\), we say that \(D\) is an *inner* derivation:

\[ D(d_1d_2) = d_1d_2f-fd_1d_2 = d_1d_2f-d_1fd_2 + d_1fd_2-fd_1d_2 = d_1(d_2f-fd_2) + (d_1f-fd_1)d_2 = d_1D(d_2) + D(d_1)d_2 \]

Dirac showed that all derivations are inner derivations for some \(f\in{\mathcal A}\). The package supports derivations:

```
f <- rweyl()
D <- as.der(f) # D(x) = xf-fx
```

Then

```
d1 <- rweyl()
d2 <- rweyl()
D(d1*d2) == d1*D(d2) + D(d1)*d2
```

`## [1] TRUE`

In the package, the product is customisable. In general, product `a*b`

[where `a`

and `b`

are `weyl`

objects] is dispatched to the following sequence of functions:

`weyl_prod_multivariate_nrow_allcolumns()`

`weyl_prod_multivariate_onerow_allcolumns()`

`weyl_prod_multivariate_onerow_singlecolumn()`

`weyl_prod_univariate_onerow()`

`weyl_prod_helper3()`

[default]

In the above, “univariate” means “generated by \(\left\lbrace x,\partial_x\right\rbrace\)” [so the corresponding `spray`

object has *two* columns]; and “multivariate” means that the algebra is generated by more than one variable, typically something like \(\left\lbrace x,y,z,\partial_x,\partial_y,\partial_z\right\rbrace\).

The penultimate function `weyl_prod_univariate_onerow()`

is sensitive to option `prodfunc`

which specifies the recurrence relation used. This defaults to `weyl_prod_helper3()`

:

`weyl_prod_helper3`

```
## function (a, b, c, d)
## {
## f <- function(r) {
## factorial(r) * choose(b, r) * choose(c, r)
## }
## ind <- numeric(0)
## val <- numeric(0)
## for (r in 0:b) {
## ind <- rbind(ind, c(a + c - r, b + d - r))
## val <- c(val, f(r))
## }
## spray(ind, val, addrepeats = TRUE)
## }
## <bytecode: 0x561fc10f09a8>
## <environment: namespace:weyl>
```

Function `weyl_prod_helper3()`

follows Wolf. This gives the univariate concatenation product \((\partial^a x^b)(\partial^c x^d)\) in terms of standard generators:

\[ \partial^a x^b \partial^c x^d=\sum_{r=0}^b r!{b\choose r}{c\choose r} \partial^{a+c-r}x^{b+d-r} \]

The package also includes lower-level function `weyl_prod_helper1()`

implementing \(\partial^a x^b \partial^c x^d=\partial^ax^{b-1}\partial^cx^{d+1} + c\partial^ax^{b-1}\partial^{c-1}x^d\) (together with suitable bottoming-out). I expected function `weyl_prod_helper3()`

to be much faster than `weyl_prod_helper1()`

but there doesn’t seem to be much difference between the two.

We can exploit this package customisability by considering, instead of \(\left\lbrace x,\partial\right\rbrace\), the algebra generated by \(\left\lbrace e,\partial\right\rbrace\), where \(e\) maps \(f\) to \(e^xf\): if \(f\) maps \(x\) to \(f(x)\), then \(ef\) maps \(x\) to \(e^xf(x)\). We see that \(\partial e-e\partial=e\). With this, we can prove that \(\partial^ne=e(1+\partial)^n\) and \(e^n\partial=e^n\partial+ne^n\) and, thus

\[ (e^a\partial^b)(e^c\partial^d) =e^{a+1}(1+\partial)^be^{c-1}\partial^d =e^{a}\partial^{b-1}e^{c}\partial^{d+1}+ce^{a}\partial^{b-1}e^{c}\partial^d\]

We may implement this set in package idiom as follows:

```
`weyl_e_prod` <- function(a,b,c,d){
if(c==0){return(spray(cbind(a,b+d)))}
if(b==0){return(spray(cbind(a+c,d)))}
return(
Recall(a,b-1,c,d+1) +
c*Recall(a,b-1,c,d) # cf: c*Recall(a,b-1,c-1,d)) for regular Weyl algebra
)
}
```

Then, for example, to calculate \(\partial^2e=e(1+2\partial+\partial^2)\):

```
options(prodfunc = weyl_e_prod)
options(weylvars = "e") # changes print method
d <- weyl(spray(cbind(0,1)))
e <- weyl(spray(cbind(1,0)))
d*d*e
```

```
## A member of the Weyl algebra:
## +e +2*e*d +e*d^2
```

`d^2*e`

```
## A member of the Weyl algebra:
## +e +2*e*d +e*d^2
```

By way of verification:

`d^5*e == e*(1+d)^5`

`## [1] TRUE`

which verifies that indeed \(\partial^5e=e(1+\partial)^5\). Another verification would be to cross-check with Mathematica, here working with \(\partial e\partial^2e\):

`In[1] := D[Exp[x]*D[D[Exp[x]*f[x],x],x],x]`

`Out[1] := 2E^2x f[x] + 5E^2x f'[x] + 4E^2xf''[x] + E^2x f'''[x]`

```
options(polyform = TRUE)
d*e*d^2*e
```

```
## A member of the Weyl algebra:
## +2*e^2 +5*e^2*d +4*e^2*d^2 +e^2*d^3
```

We can manipulate more complicated expressions too. Suppose we want to evaluate \((1+e^2\partial)(1-5e^3\partial^3)\):

```
o1 <- weyl(spray(cbind(2,1)))
o2 <- weyl(spray(cbind(3,3)))
options(polyform = FALSE)
(1+o1)*(1-5*o2)
```

```
## A member of the Weyl algebra:
## e d val
## 5 3 = -15
## 3 3 = -5
## 5 4 = -5
## 2 1 = 1
## 0 0 = 1
```

And of course we can display the result in symbolic form:

```
options(polyform = TRUE)
(1+o1)*(1-5*o2)
```

```
## A member of the Weyl algebra:
## 1 -15*e^5*d^3 -5*e^3*d^3 -5*e^5*d^4 +e^2*d
```

S. C. Coutinho 1997.

*The many avatars of a simple algebra*. The American Mathematical Monthly, 104(7):593-604. DOI https://doi.org/10.1080/00029890.1997.11990687.K. B. Wolf 1975.

*The Heisenberg-Weyl ring in quantum mechanics*. In*Group theory and its applications*, Volume III, Academic Press, Inc