Manipulation of sympy Functions

Question

I am working over arbitrary discrete functions like x₁(k), x₂(k), etc. where k is an Integer, and for any k, x_i(k) is a real number (however, I might add that for the symbolic computations I want to perform, this is not relevant).

I need to manipulate those functions symbolically, i.e., I need to be able to perform typical operations like multiplication or addition of these functions (e.g. x₁(k)x₂(k)) as well as manipulation of their arguments (e.g. x₁(2k+1)). Due to the argument manipulation, I can not just use sympy's Symbol class because they do not have arguments. Instead, sympy's Function class supports arguments and I thought that is exactly what I need. But for things like x1*x2 it throws

TypeError: unsupported operand type(s) for *: 'UndefinedFunction' and 'UndefinedFunction'

Thus, I am not quite sure if the Function class is even suited for what I am looking for, or if I am just misunderstanding how to achieve what I need. How can I use sympy for these kind of computations?

EDIT: For clarification, (before I knew that the cls option existed in sympy.Symbol) I thought that

import sympy as sp

x1 = sp.Symbol("x1")

wasn't going to work, due to the lack of arguments. So I tried

k = sp.Symbol("k")
x1 = sp.Function("x1")(k)

which would have arguments, but throw the error for things like x1*x1. Hope that is more clear.

Follow up question: As suggested in the answer of ForceBru, if I use

k = sp.Symbol("k")
x1 = sp.symbols("x1", cls=sp.Function)(k)
x2 = sp.symbols("x2", cls=sp.Function)(k)

then this seems to do the trick. However, what I also need to do is raise the argument by an integer j, i.e., I might have an expression like 2x₁(2k+1)x₂(2k+3) (or even more complicated) which after applying this operation returns 2x₁(2k+1+j)x₂(2k+3+j). Can't really use expr.args since that would not only return the arguments of x₁ and x₂. Also, expr.subs may likely be risky to use if not only arguments of x₁ and x₂ are affected.

Do you see any way how this can be achieved?

EDIT2: After playing around more I have come to what seems to work by using expr.atoms(sp.Function) which returns a list of functions which can then be replaced by (j-times) shifted versions of themselves. Will accept FroceBru's answer as this was the crucial hint.

Answer 1

Your code tries to multiply the function themselves, not the result of their application.

I need to be able to perform typical operations like multiplication or addition of these functions (e.g. x1(k)x2(k))

You can do this like this:

>>> import sympy as sp
>>> x1, x2 = sp.symbols("x_1 x_2", cls=sp.Function)
>>> k = sp.Symbol('k')
>>> x1(k) * x2(k) # simple multiplication
x_1(k)*x_2(k)
>>> x1(2*k + 1) * x2(3*k + 2) # argument manipulation
x_1(2*k + 1)*x_2(3*k + 2)

The difference is that x1 * x2 attempts to multiply the functions themselves, but it doesn't make much sense, unless * means convolution. For example, what's sin * cos? On the contrary, x1(k) * x2(k) is a reasonable expression, just like sin(x) * cos(x), because it's a product of the results of applying x1 and x2 to k.