Summing vector products in SymPy

Question

The following code computes $\sum_{i=1}^3 x_i y_i$ symbolically first, and then substitutes specific numeric values into the expression to yield a single number. (I've added some lines to show intermediate results as well.)

IPython console for SymPy 1.4 (Python 3.7.4-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

In [1]: i = symbols('i')                                                                                                
In [2]: x = IndexedBase('x')                                                                                            
In [3]: y = IndexedBase('y')                                                                                            
In [4]: xes = [(x[1], 9), (x[2], 8), (x[3], 7)]                                                                         
In [5]: ys = [(y[1], 4), (y[2], 5), (y[3], 6)]                                                                          
In [6]: Snum = Sum(x[i] * y[i],(i,1,3))                                                                                 
In [7]: Snum.doit()                                                                                                     
Out[7]: x[1]⋅y[1] + x[2]⋅y[2] + x[3]⋅y[3]
In [8]: Sx = Snum.doit().subs(xes)                                                                                      
In [9]: Sx                                                                                                              
Out[9]: 9⋅y[1] + 8⋅y[2] + 7⋅y[3]
In [10]: Sxy = Sx.subs(ys)                                                                                              
In [11]: Sxy                                                                                                            
Out[11]: 118

I've been trying to transfer this idea to a case which involves vectors instead of scalars and cross products instead of multiplication, but without success.

In [12]: from sympy.vector import *                                                                                     
In [13]: N = CoordSys3D('N')                                                                                            
In [14]: r = IndexedBase('r')                                                                                           
In [15]: F = IndexedBase('F')                                                                                           
In [16]: rs = [(r[1], N.i + 2*N.j + 4*N.k), (r[2], 3*N.i - 3*N.j + 4*N.k), (r[3], -N.i + 5*N.j + 2*N.k)]                
In [17]: Fs = [(F[1], -2*N.i - 1*N.j + 4*N.k), (F[2], -3*N.i - 2*N.j + 1*N.k), (F[3], N.i - 2*N.j - 3*N.k)]             
In [18]: rs[0][1] ^ Fs[0][1]                                                                                            
Out[18]: 12*N.i + (-12)*N.j + 3*N.k
In [19]: Svec = Sum(r[i] ^ F[i],(i,1,3))                                                                                
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-19-08e9d4fdd092> in <module>
----> 1 Svec = Sum(r[i] ^ F[i],(i,1,3))

TypeError: unsupported operand type(s) for ^: 'Indexed' and 'Indexed'

It seems indexed objects do not support vector operations like cross product.

My question is, is there any way in SymPy to perform a computation similar to my purely numeric example above which involves vectors and vector operations? I'd like to compute $T = \sum_{i=1}^3 r_i x F_i$ in such a way that first I get r_1 x F_1 + r_2 x F_2 + r_3 x F_3, and then I can make the appropriate substitutions to get a single vector 6*N.i + (-28)*N.j + (-15)*N.k as a result. Is this possible? If so, how?

Answer 1

Can you do the sum as a product but then result.replace(lambda x: x.is_Mul, lambda x: reps[x.args[0]]^reps[x.args[1]]) where reps = dict(rs); reps.update(dict(Fs))? e.g. using ints and bitwise operations works:

>>> rs = [(r[1], 1), (r[2], 2), (r[3], 3)]
>>> Fs = [(F[1], 4), (F[2], 5), (F[3], 6)]
>>> Sum(r[i]*F[i], (i, 1, 3)).doit().replace(
...    lambda x: x.is_Mul,
...    lambda x: int(reps[x.args[0]])^int(reps[x.args[1]]))
... 
17

A problem with this, however, is that (currently) the IndexedBase cannot be declared as noncommutative. Instead of IndexedBase you can use Function:

>>> f=Function('f', commutative=False)
>>> g=Function('g', commutative=False)
>>> Sum(f(i)*g(i), (i, 1, 3)).doit()
f(1)*g(1) + f(2)*g(2) + f(3)*g(3)
>>> Sum(g(i)*f(i), (i, 1, 3)).doit()
g(1)*f(1) + g(2)*f(2) + g(3)*f(3)

In your case you will name them as r and F instead of f and g and will map *them( to the vectors like rs = [(r(1), N.i + 2*N.j + 4*N.k), ....

If the following diff were applied to SymPy, then I think the IndexedBase declared as commutative=False would work:

diff --git a/sympy/tensor/indexed.py b/sympy/tensor/indexed.py
index 9024f21..f02aa3e 100644
--- a/sympy/tensor/indexed.py
+++ b/sympy/tensor/indexed.py
@@ -138,7 +138,6 @@ class Indexed(Expr):
     True

     """
-    is_commutative = True
     is_Indexed = True
     is_symbol = True
     is_Atom = True
@@ -419,7 +418,6 @@ class IndexedBase(Expr, NotIterable):
     >>> C_inherit == C_explicit
     True
     """
-    is_commutative = True
     is_symbol = True
     is_Atom = True

For example:

>>> I=IndexedBase('i', commutative=False)
>>> I[1]*I[2]
i[1]*i[2]
>>> I[2]*I[1]
i[2]*i[1]

You might want to open an issue for this here.