Problemi durante la sostituzione di una matrice in un polinomio

Question

Esempio: let

M = Matrix([[1,2],[3,4]]) # and 
p = Poly(x**3 + x + 1)    # then
p.subs(x,M).expand()

dà l'errore:

TypeError: impossibile aggiungere <class'sympy.matrices.immutable.ImmutableDenseMatrix '> e <class' sympy.core.numbers.One '>

il che è molto plausibile poiché i due primi termini diventano matrici ma l'ultimo termine (il termine costante) non è una matrice ma uno scalare. Per rimediare a questa situazione ho cambiato il polinomio in

p = Poly(x**3 + x + x**0)    # then

lo stesso errore persiste. Sono obbligato a digitare l'espressione a mano, sostituendo x con M? In questo esempio il polinomio ha solo tre termini ma in realtà incontro (polinomi multivariati con) dozzine di termini.

Answer 1

Quindi penso che la domanda ruoti principalmente attorno al concetto di polinomio Matrix :

(dove P è un polinomio e A è una matrice)

Penso che questo significhi che il termine libero è un numero e non può essere aggiunto con il resto che è una matrice, in effetti l'operazione di addizione è indefinita tra questi due tipi.

TypeError: impossibile aggiungere <class'sympy.matrices.immutable.ImmutableDenseMatrix '> e <class' sympy.core.numbers.One '>

Tuttavia, questo può essere aggirato definendo una funzione che valuta il polinomio della matrice per una matrice specifica. La differenza qui è che stiamo usando l' esponenziazione della matrice , quindi calcoliamo correttamente il termine libero del polinomio della matrice a_0 * Idove si I=A^0trova la matrice identità della forma richiesta:

from sympy import *
x = symbols('x')
M = Matrix([[1,2],[3,4]])
p = Poly(x**3 + x + 1)

def eval_poly_matrix(P,A):
    res = zeros(*A.shape)
    for t in enumerate(P.all_coeffs()[::-1]):
        i, a_i = t
        res += a_i * (A**i)
    return res

eval_poly_matrix(p,M)

Produzione:

In questo esempio il polinomio ha solo tre termini ma in realtà incontro (polinomi multivariati con) dozzine di termini.

La funzione eval_poly_matrixsopra può essere estesa per funzionare per polinomi multivariati utilizzando il .monoms()metodo per estrarre monomi con coefficienti diversi da zero , in questo modo:

from sympy import *
x,y = symbols('x y')
M = Matrix([[1,2],[3,4]])
p = poly( x**3 * y + x * y**2 + y )

def eval_poly_matrix(P,*M):
    res = zeros(*M[0].shape)
    for m in P.monoms():
        term = eye(*M[0].shape)
        for j in enumerate(m):
            i,e = j
            term *= M[i]**e
        res += term
    return res

eval_poly_matrix(p,M,eye(M.rows))

Nota: sono possibili alcuni controlli di integrità, gestione dei casi limite e ottimizzazioni:

1. Il numero di variabili presenti nel polinomio si riferisce al numero di matrici passate come parametri (la prima non dovrebbe mai essere maggiore della seconda, e se è inferiore a qualche logica deve essere presente per gestirla, ho gestito solo il caso quando i due sono uguali)
2. Tutte le matrici devono essere quadrate secondo la definizione del polinomio di matrice
3. Una discussione su una versione multivariata della regola di Horner compare nei commenti di questa domanda . Ciò potrebbe essere utile per ridurre al minimo il numero di moltiplicazioni di matrici.
4. Gestire il fatto che in una matrice il polinomio x*yè diverso dal y*xperché la moltiplicazione della matrice non è commutativa . Apparentemente le funzioni poly in sympy non supportano variabili non commutative , ma puoi definire simboli con commutative=Falsee sembra esserci una via da seguire

Circa il 4 ° punto sopra, c'è il supporto per le espressioni Matrix in SymPy e questo può aiutare qui:

from sympy import *
from sympy.matrices import MatrixSymbol

A = Matrix([[1,2],[3,4]])
B = Matrix([[2,3],[3,4]])

X = MatrixSymbol('X',2,2)
Y = MatrixSymbol('Y',2,2)
I = eye(X.rows)

p = X**2 * Y + Y * X ** 2 + X ** 3 - I
display(p)

p = p.subs({X: A, Y: B}).doit()
display(p)

Produzione:

Per ulteriori sviluppi su questa funzione, segui # 18555