Подставляем матрицу Sympy в полином

У меня есть матрица Sympy A и полиномиальное выражение P, и я хотел бы вычислить P (A).

Вот пример:

x = Symbol('x')
P = x**2 - 3*x + 5
A = Matrix([ [1,3], [-1,2] ])
P.subs(x, A)

Я ожидаю, что Sympy вычислит A**2 - 3*A + 5*eye(2) (в этом примере результатом является нулевая матрица).

Но это не срабатывает с сообщением об ошибке:

AttributeError: ImmutableMatrix has no attribute as_coeff_Mul.

Есть ли способ получить то, что я хочу?

Изменить: Я попытался преобразовать P в класс полиномов Sympy, а затем заменить его, но результат бесполезен:

Poly(P).subs(A)
Poly(Matrix([ [ 1, 3], [-1, 2]])**2 - 3*Matrix([ [ 1, 3],        
    [-1, 2]]) + 5, Matrix([ [ 1, 3], [-1, 2]]), domain='ZZ')

Я могу получить правильный результат с помощью следующей функции:

def poly_matrix(P, A):
    coeffs = Poly(P).all_coeffs()[::-1]
    res = zeros(A.rows)
    for i in range(len(coeffs)):
        res += coeffs[i]*(A**i)
    return res

Но все же ищу более производительный встроенный вариант.


python sympy matrix polynomials
person Adrien    schedule 07.09.2015    source источник

Ответы (3)


arrow_upward
2
arrow_downward

Матричные выражения все еще содержат некоторые ошибки, надеюсь, они будут исправлены в будущем.

В любом случае, вот альтернативный способ выполнить замену:

In [1]: x = MatrixSymbol('x', 2, 2)

In [2]: P = x**2 - 3*x + 5*eye(2)

In [3]: P
Out[3]: 
                 2
⎡5  0⎤ + -3⋅x + x 
⎢    ⎥            
⎣0  5⎦            

In [4]: A = Matrix([ [1,3], [-1,2] ])

In [5]: P.subs(x, A)
Out[5]: 
                               2
⎡5  0⎤ + -3⋅⎡1   3⎤ + ⎛⎡1   3⎤⎞ 
⎢    ⎥      ⎢     ⎥   ⎜⎢     ⎥⎟ 
⎣0  5⎦      ⎣-1  2⎦   ⎝⎣-1  2⎦⎠ 

In [6]: P.subs(x, A).doit()
Out[6]: 
                   2
⎡2  -9⎤ + ⎛⎡1   3⎤⎞ 
⎢     ⎥   ⎜⎢     ⎥⎟ 
⎣3  -1⎦   ⎝⎣-1  2⎦⎠

Здесь оказывается, что MatPow не может выполнить операцию .doit (). Наверное, это очередной баг.

В выводе номер 5 также есть ошибка с принтером (+ -3 должно быть просто -3).

Очень надеюсь, что кто-то со временем исправит все эти проблемы.

person Francesco Bonazzi    schedule 15.09.2015
comment
Интересная альтернатива. Ошибка doit() исправлена ​​в 0.7.7.dev. - person Adrien; 18.09.2015

arrow_upward
0
arrow_downward

Попробуйте вычислить каждый член многочлена.

(x**2).subs(x,A) - (3*x).subs(x,A) + 5*(eye(2))

Это оценит ваше выражение.

person rresol    schedule 17.09.2015
comment
Что ж, если бы я заменял термины один за другим вручную, я бы также с самого начала набрал A**2 - 3*A + 5*eye(2). - person Adrien; 18.09.2015
comment
буду искать лучшее решение. - person rresol; 19.09.2015

arrow_upward
0
arrow_downward

Чтобы оценить P на A, вы можете заменить x на Matrix и преобразовать постоянный член умножением на eye(2):

P_ = Poly(P, x)
(P_ - P_.coeff_monomial(1)).as_expr().subs(x, A) * eye(2) + P_.coeff_monomial(1) * eye(2) # P(A)

Первое умножение на eye(2) гарантирует, что первый член в сумме является матрицей, даже если P является просто константой, то есть P_ - P_.coeff_monomial(1) == 0.

person t.y    schedule 12.07.2020