How to generate a symbolic multivariate polynomial of a given dimension in SymPy?

603 views
Skip to first unread message

foadsf

unread,
Aug 1, 2018, 7:01:13 AM8/1/18
to sympy

I want to use power series to approximate some PDEs. The first step I need to generate symbolic multivariate polynomials, given a numpy ndarray.

Consider the polynomial below:


I want to take a m dimensional ndarray of D=[d1,...,dm] where djs are non-negative integers, and generate a symbolic multivariate polynomial in the form of symbolic expression. The symbolic expression consists of monomials of the form:


Fo example if D=[2,3] the output should be


For this specific case I could nest two for loops and add the expressions. But I don't know what to do for Ds with arbitrary length. If I could generate the D dimensional ndarrays of A and X without using for loops, then I could use np.sum(np.multiply(A,X)) as Frobenius inner product to get what I need.

I would appreciate if you could help me know how to do this in SymPy. Thanks in advance.

Kalevi Suominen

unread,
Aug 1, 2018, 9:10:35 AM8/1/18
to sympy
I would first generate the list of monomial indices by using e.g. itertool.product, then create a dictionary containing the indexed coefficients, and finally create a Poly object with given variables from those coefficients. For example, to construct a polynomial of degree 3 or less in 3 variables this could be done:

indices = [i for i in itertools.product(range(4), repeat=3) if sum(i) < 4]
a = IndexBase('a')
coeffs = {i: a[i] for i in indices}
vars = symbols('x:3')
Poly(coeffs, *vars)

Kalevi Suominen

Kalevi Suominen

unread,
Aug 1, 2018, 9:12:20 AM8/1/18
to sympy
Correction: IndexedBase instead of IndexBase

Foad Sojoodi Farimani

unread,
Aug 1, 2018, 9:42:20 AM8/1/18
to sy...@googlegroups.com
Dear Kalevi,

First of thanks a lot for the great instructions. It is a good step forward. I applied your code as:

    from itertools import product
    from sympy import IndexedBase, symbols, Poly
    d=3
    l=4
    indices = [i for i in product(range(l), repeat=d) if sum(i) < l]
    a = IndexedBase('a')
    coeffs = {i: a[i] for i in indices}
    vars = symbols('x:3')
    Poly(coeffs, *vars)


However there are a couple of issues:
  • `  if sum(i) < l ` causes the total degree of each monomial to be less than 4, that's not what I want. I want to take a list or ndarray of non-negative integers D=[d1,...,dm] and each monomial should have a degree of  i_j for x_j which is less than d_j , so not the total degree of each monomial but the degree of each variable is important. is there any way to tell `itertools.product` to give the Cartesian product of a list of ranges. something like  itertools.product(range(d1),...,range(dm)). 
  • in the second line from bottom how can I change the `3` with a variable?
Thanks for your help again.

Best,
Foad




--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/1c70ac48-9c09-43b5-9374-15acb37d2860%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Foad Sojoodi Farimani

unread,
Aug 1, 2018, 10:20:26 AM8/1/18
to sy...@googlegroups.com
I think I solved the second problem. From here I learned that if I change `symbols('x:3')` to `symbols(f'x1:{d+1}')` the second issue is solved.

Foad Sojoodi Farimani

unread,
Aug 1, 2018, 11:28:43 AM8/1/18
to sy...@googlegroups.com
From here I could solve the problem:

    from itertools import product
    from sympy import IndexedBase, symbols, Poly
    D=[2,3]
    d=len(D)
    indices=list(product(*map(range, D)))
    a = IndexedBase('a')
    coeffs = {i: a[i] for i in indices}
    vars = symbols(f'x1:{d+1}')
    Poly(coeffs, *vars)

and the result is:

    Poly(a[1, 2]*x1*x2**2 + a[1, 1]*x1*x2 + a[1, 0]*x1 + a[0, 2]*x2**2 + a[0, 1]*x2 + a[0, 0], x1, x2, domain='EX')

Foad Sojoodi Farimani

unread,
Aug 1, 2018, 11:43:07 AM8/1/18
to sy...@googlegroups.com
other solutions have been posted here and here 
Reply all
Reply to author
Forward
0 new messages