How to find Coefficient of more than one variable of desired power in Python

Question

if I had an equation with more than one variable, lets say x**2+x*y+y**2, I could find Coefficients with CoefficientList from Mathematica. It would print the desired matrix as 3x3.

In Python, I found sympy.coeff(x, n) but I can not return the coefficient for more than one variable.

Is there a way you know?

Does this answer help? https://stackoverflow.com/questions/22955888/how-to-extract-all-coefficients-in-sympy — Marcus, Mar 04 '20 at 12:11
Does this answer your question? [How to extract all coefficients in sympy](https://stackoverflow.com/questions/22955888/how-to-extract-all-coefficients-in-sympy) — KJTHoward, Mar 04 '20 at 12:12
@Marcus This post is for coefficient of one variable polys only. IF you look [here](https://reference.wolfram.com/language/ref/CoefficientList.html?fbclid=IwAR1BQnChPQLsjQKebxjnwDIWR9mgUXKyBSp1NAn0mZNzTWNyMVI93zokbTM) the 3rd example you will see what i would love to have as output (its more than one variable) — Billy Matlock, Mar 04 '20 at 12:24
@KJTHoward Didn't let my quote you to, so the previous answer holds, as its the same post :) — Billy Matlock, Mar 04 '20 at 12:25
It seem you can use the monomial methods from sympy. Have a look at this answer. https://stackoverflow.com/questions/36818025/best-way-to-isolate-one-coefficient-of-a-multivariate-polynomial-in-sympy — Marcus, Mar 04 '20 at 12:47
Really thank you! Don't know why I didn't find this answer before. Gonna take a look at it and the Answer below. Hope i'll understand them :) — Billy Matlock, Mar 04 '20 at 16:59

score 2 · Answer 1 · answered Mar 04 '20 at 17:25

For the quadratic case, the following matches the output of Mathematica:

def CoefficientList(p, v):
    """
    >>> CoefficientList(x**2+2*x*y+3*y**2+4*x+5*y+6,(x,y))
    Matrix([
    [6, 5, 3],
    [4, 2, 0],
    [1, 0, 0]])
    """
    assert len(v) == 2
    n = len(v)
    t = [prod(i) for i in subsets([S.One] + list(v), n, repetition=n)]
    p = p.expand()
    m = zeros(n + 1)
    r = n + 1
    while t:
        if r > n:
            n -= 1
            r = 0
            c = n
        x = t.pop()
        if x == 1:
            d = p
        else:
            p, d = p.as_independent(x, as_Add=True)
        co = d.as_coeff_Mul()[0]
        m[r, c] = co
        r += 1
        c -= 1
    return m

But the monomials method is a good one to use. To get all the coefficients of all monomials I would recommend storing them in a defaultdict with a default value of 0. You can then retrieve the coefficients as you wish:

def coeflist(p, v):
    """
    >>> coeflist(x**2+2*x*y+3*y**2+4*x+5*y+6, [x])
    defaultdict(<class 'int'>, {x**2: 1, x: 2*y + 4, 1: 3*y**2 + 5*y + 6})
    >>> coeflist(x**2+2*x*y+3*y**2+4*x+5*y+6, [x, y])
    defaultdict(<class 'int'>, {x**2: 1, x*y: 2, x: 4, y**2: 3, y: 5, 1: 6})
    >>> _[y**2]
    3
    """
    from collections import defaultdict
    p = Poly(p, *v)
    rv = defaultdict(int)
    for i in p.monoms():
        rv[prod(i**j for i,j in zip(p.gens, i))] = p.coeff_monomial(i)
    return rv

JohanC · Accepted Answer · 2020-03-04T14:22:47.650

Here is a way to find each of the coefficients of the quadratic form and represent them as a matrix:

import sympy as sy
from sympy.abc import x, y

def quadratic_form_matrix(expr, x, y):
    a00, axx, ayy, ax, ay, axy = sy.symbols('a00 axx ayy ax ay axy')
    quad_coeffs = sy.solve([sy.Eq(expr.subs({x: 0, y: 0}), a00),
                            sy.Eq(expr.diff(x).subs({x: 0, y: 0}), 2 * ax),
                            sy.Eq(expr.diff(x, 2).subs({y: 0}), 2 * axx),
                            sy.Eq(expr.diff(y).subs({x: 0, y: 0}), 2 * ay),
                            sy.Eq(expr.diff(y, 2).subs({x: 0}), 2 * ayy),
                            sy.Eq(expr.diff(x).diff(y), 2 * axy)],
                           (a00, axx, ayy, ax, ay, axy))
    return sy.Matrix([[axx, axy, ax], [axy, ayy, ay], [ax, ay, a00]]).subs(quad_coeffs)

expr = x**2 + 2*x*y + 3*y**2 + 7
M = quadratic_form_matrix(expr, x, y)
print(M)

XY1 = sy.Matrix([x, y, 1])
quadatric_form = (XY1.T * M * XY1)[0]
print(quadatric_form.expand())

PS: Applying a conversion to a multivariate polygon as suggested @Marcus' reference, and then converting to a matrix would result in following code. Note that to get the constant term, 1 can be passed to coeff_monomial(1). The non-diagonal elements of the symmetric matrix for the quadratic form need to be half of the corresponding coefficients.

import sympy as sy
from sympy.abc import x, y

def quadratic_form_matrix_using_poly(expr, x, y):
    p = sy.poly(expr)
    axx = p.coeff_monomial(x * x)
    ayy = p.coeff_monomial(y * y)
    a00 = p.coeff_monomial(1)
    ax = p.coeff_monomial(x) / 2
    ay = p.coeff_monomial(y) / 2
    axy = p.coeff_monomial(x * y) / 2
    return sy.Matrix([[axx, axy, ax], [axy, ayy, ay], [ax, ay, a00]])

expr = x**2 + 2*x*y + 3*y**2 + 7 + 11*x + 23*y

M = quadratic_form_matrix(expr, x, y)
print(M)

XY1 = sy.Matrix([x, y, 1])
quadatric_form = (XY1.T * M * XY1)[0]
print(quadatric_form.expand())

really appreciate it! Gonna study it and adapt it to my system (its a bit more complex due to Constants and Functions). If I have a problem, may I add a comment to this Answer or nah? — Billy Matlock, Mar 04 '20 at 17:00

How to find Coefficient of more than one variable of desired power in Python

2 Answers2