Hi,
I am trying to write a symbolic library for computing 3D polynomials in an analytical fashion (the variable is a real value t, and the monomials of the polynomials are 3D vectors). In particular, I want to compute the cross product of two polynomials (in a follow up to that question ). In Sympy, I found that I can either:
Is there any way that I can represent symbolically the cross product in sympy?
edit: In particular, I am interested in nullifying the cross product of two identical vectors, and taking into account the anticommutative property when factorizing a polynomial in terms of one of the vectors.
new edit: To make myself more clear, I want to stay on a "symbolic level". That is, I don't want to develop the terms of my vector along each variable.
For instance, here is the code to compute a Bezier curve:
from sympy import *
init_printing(use_unicode=True)
from scipy.special import binom
#default control points of a Bezier curve
p_is = [symbols('p'+str(i)) for i in range(5)]
class Bezier:
#eq is the equation of the curve, pis are the stored control point
#for special purpose such as the derivatives
def __init__(self, deg, eq = None, pis = None):
assert(deg)
self.deg = deg;
n_pts = deg +1
if pis == None:
self.pis = p_is[:n_pts]
else:
self.pis = pis[:];
if eq == None:
#computing the bernstein polynoms for a given degree
factors = [binom(deg,i) * t**i * (1-t)**(deg-i)*pis[i] for i in range(n_pts)]
self.eq = sum(factors);
else:
self.eq = eq
def __repr__(self):
res = "Degree : " + str(self.deg)
res += "\nEquation : " + self.eq.__repr__()
res += "\nwaypoints :\n" + str(self.pis)
return res
def __str__(self):
return self.__repr__()
b = Bezier(3)
print b
# Degree : 3
# Equation : 1.0*p0*(-t + 1)**3 + 3.0*p1*t*(-t + 1)**2 + 3.0*p2*t**2*(-t + 1) + 1.0*p3*t**3
# waypoints :
# [p0, p1, p2, p3]
print b.eq
# 1.0*p0*(-t + 1)**3 + 3.0*p1*t*(-t + 1)**2 + 3.0*p2*t**2*(-t + 1) + 1.0*p3*t**3
As one can see, the fact that the variables p_is are vector is not really relevant, except in the case where the cross product appears. If I were to compute the cross product of b with itself, then several terms would disappear because the some vectors would be crossed against themselves.
What I tried to do was to "emulate" the cross product with a simple multiplication, and then iterate over the resulting equation to remove all the squared terms (that correspond to zero). But this is not sufficient because the product does not preserve the anticomutative aspect of the cross product.
What I would really like would be for the actual cross product sign (say X) appear in the equation. I hope I am clearer
Thanks a lot for your help
Steve
I found a solution that works for my case (where I only apply cross product once, some reflexion is needed to consider an extension). The idea is to use additional symbols to represent the skew symmetric-matrix that achieves the cross product. For instance p0 ^ p1 will be written Cp0 * p1. I implemented a method that achieves this, and also makes sure that if pi ^ pi = 0.
To allow a better factorization, I introduced an arbitrary order of transformation based on the alphabetic order of the symbols. This means that p2 ^ p1 will be written -Cp1 * p2 The reason for this is that otherwise simplifications such as Cp2 ^ p1 + Cp1 ^p2 = 0 will not be detected by sympy.
Anyway, this is rather hacky but in my case it allowed me to write my library, so here it is.
The method that performs the cross product is "cross", at the end of the file.
#declaring vector symbols variables
p_is = [symbols('p'+str(i)) for i in range(20)]
#the cross product will be represented
#by the skew symmetric matrix Cpi for a vector pi.
#Since in the resulting equation, symbols seem to appeer
#in alphabetic order, the matrix multiplication will be coherent
p_isX = [symbols('Cp'+str(i)+'') for i in range(20)]
#achieves the cross product between two elements
#s0 and s1 are the considered vector symbols (here, two pis)
#e0 and e1 are the scalar multiplier of each term
def crossElement(s0,e0,s1,e1):
if s0 == s1:
return 0
else:
# here, take order of apparition into account to allow
# future factorization. Otherwise
# something like p0 X p1 + p1 X p0 will not be dientified as zero
id0 = p_is.index(s0)
id1 = p_is.index(s1)
if(id0 < id1):
return simplify(e1*e0*p_isX[id0]*s1)
else:
return simplify((-e1)*e0*p_isX[id1]*s0)
#retrieve all the pis and associate scalar factors
#from an equation
def getAllSymbols(eq):
res = []
for el in p_is:
dic = eq.as_poly(el).as_dict()
if(len(dic) > 1):
res += [(el, dic[(1,)])]
return res;
#generates the cross product of two equations,
#where each pi term is a vector, and others
#are scalar variable.
#cross product
def cross(eq1,eq2):
decomp1 = getAllSymbols(eq1)
decomp2 = getAllSymbols(eq2)
res = 0
#applies distributive cross product between
#all terms of both equations
for (s0, e0) in decomp1:
for (s1, e1) in decomp2:
res += crossElement(s0,e0,s1,e1)
return res