Vector expressions with SymPy#
First, install and import Vector and SymPy.
[1]:
import vector
import sympy
How the SymPy backend differs from the others#
SymPy is a computer algebra system like Mathematica and Maple. It primarily deals with algebraic expressions, rather than concrete numbers. However, all of the coordinate transformations and vector manipulations can be applied symbolically through Vector’s SymPy backend.
When comparing SymPy to the other backends, note that SymPy vector expressions have a different sign convention for operations on space-like and negative time-like 4D vectors. For all other backends, Vector’s conventions were chosen to agree with popular HEP libraries, particularly ROOT, but for the SymPy backend, those conventions would insert piecewise if-then branches, which would complicate symbolic expressions.
When vector expressions are evaluated numerically, you can expect agreement in 2D and 3D vector operations, as well as 4D vector operations if all of the vectors have a positive time-like part (which is necessary for real momentum vectors and causal relationships between events).
Making a vector expression#
Before making a vector expression, we need symbols for each of the components, so use the sympy.symbols function. Be sure to tell SymPy to assume that they are all real-valued, not complex numbers.
[2]:
x, y, z, t, px, py, pz, eta, tau = sympy.symbols("x y z t px py pz eta tau", real=True)
Now we can make vectors just as we did with objects, though these lack concrete numerical values.
[3]:
vector.VectorSympy2D(x=x, y=y)
[3]:
VectorSympy2D(x=x, y=y)
[4]:
vector.MomentumSympy3D(px=px, py=py, pz=pz)
[4]:
MomentumSympy3D(px=px, py=py, pz=pz)
[5]:
vector.VectorSympy4D(x=x, y=y, eta=eta, tau=tau)
[5]:
VectorSympy4D(x=x, y=y, eta=eta, tau=tau)
Using a vector expression#
All of the vector operations performed on these expressions return symbolic results.
[6]:
v = vector.VectorSympy2D(x=x, y=y)
[7]:
v.rho
[7]:
[8]:
sympy.Eq(v.rho, abs(v))
[8]:
[9]:
v = vector.VectorSympy4D(x=x, y=y, z=z, t=t)
[10]:
v.tau
[10]:
[11]:
v.is_timelike()
[11]:
[12]:
boosted = v.boost(v.to_beta3())
boosted
[12]:
VectorSympy4D(x=x*(1 + x**2/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))) + x/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2) + x*y**2/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)) + x*z**2/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)), y=y*(1 + y**2/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))) + y/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2) + x**2*y/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)) + y*z**2/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)), z=z*(1 + z**2/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))) + z/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2) + x**2*z/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)) + y**2*z/(t**2*(1 + 1/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2))*(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)), t=t/sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2) + x**2/(t*sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)) + y**2/(t*sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)) + z**2/(t*sqrt(1 - x**2/t**2 - y**2/t**2 - z**2/t**2)))
[13]:
boosted.t
[13]:
They can be simplified:
[14]:
boosted.t.simplify()
[14]:
And the symbols can be replaced with numerical values:
[15]:
boosted.t.subs({x: 3, y: 2, z: 1, t: 10})
[15]:
Or converted into code for a programming language:
[16]:
import sympy.printing.fortran
print(sympy.printing.fortran.fcode(boosted.t.simplify()))
t*sqrt((t**2 - x**2 - y**2 - z**2)/t**2)*(t**2 + x**2 + y**2 + z**
@ 2)/(t**2 - x**2 - y**2 - z**2)