Identità di Newton in Python: Calcolare Somme di Potenze con SymPy senza Risolvere il Polinomio

Quando compare un polinomio di grado elevato, il riflesso naturale è quasi sempre lo stesso: trovare le radici.

Ma spesso è la strada più scomoda.

Esistono situazioni in cui non ci interessa conoscere le radici una per una.
Magari vogliamo soltanto calcolare quantità come:

[math]
r_1^5+r_2^5+r_3^5
[/math]

oppure:

[math]
r_1^{20}+r_2^{20}+r_3^{20}
[/math]

Farlo direttamente può diventare rapidamente ingestibile.

Ed è qui che entrano in gioco le identità di Newton: relazioni algebriche che collegano le somme di potenze delle radici ai coefficienti del polinomio.

La parte interessante è che queste formule si prestano perfettamente all’automazione.

Con poche righe di Python e la libreria SymPy possiamo costruire ricorrenze simboliche, generare formule automaticamente e verificare numericamente i risultati senza risolvere esplicitamente il polinomio.

In questo tutorial vedremo:

• come funzionano le identità di Newton
• come implementarle in Python
• come automatizzare le ricorrenze
• come usare il calcolo simbolico
• come verificare numericamente i risultati
• perché questa tecnica è importante nell’algebra computazionale

Tutto il codice è compatibile con Google Colab.

import sympy as sp

e1, e2, e3 = sp.symbols('e1 e2 e3')

S1, S2, S3 = sp.symbols('S1 S2 S3')

S4 = e1*S3 - e2*S2 + e3*S1

print(S4)

S1*e3 - S2*e2 + S3*e1

def newton_sums(e1, e2, e3, S1, S2, S3, n_max):
    
    sums = {
        1: S1,
        2: S2,
        3: S3
    }

    for n in range(4, n_max + 1):
        
        sums[n] = (
            e1*sums[n-1]
            - e2*sums[n-2]
            + e3*sums[n-3]
        )

    return sums

sums = newton_sums(
    e1=1,
    e2=-sp.Rational(1,2),
    e3=sp.Rational(1,6),
    S1=1,
    S2=2,
    S3=3,
    n_max=8
)

for k, v in sums.items():
    print(f"S_{k} =", sp.simplify(v))

Output atteso

S_1 = 1
S_2 = 2
S_3 = 3
S_4 = 25/6
S_5 = 6
S_6 = 301/36
S_7 = 665/54
S_8 = 1705/108

x = sp.Symbol('x')

P = x**3 - 2*x**2 - x + 1

print(P)

x**3 - 2*x**2 - x + 1

Estrazione automatica dei coefficienti

poly = sp.Poly(P, x)

coeffs = poly.all_coeffs()

print(coeffs)

[1, -2, -1, 1]

e1 = 2
e2 = -1
e3 = -1

Codice

S1 = e1

S2 = e1**2 - 2*e2

S3 = e1*S2 - e2*S1 + 3*e3

print(S1, S2, S3)

e1 e1**2 - 2*e2 -e1*e2 + e1*(e1**2 - 2*e2) + 3*e3

Generazione completa automatica

sums = newton_sums(
    e1,
    e2,
    e3,
    S1,
    S2,
    S3,
    10
)

for n, val in sums.items():
    print(f"S_{n} =", sp.expand(val))

S_1 = S1
S_2 = S2
S_3 = S3
S_4 = S1*e3 - S2*e2 + S3*e1
S_5 = S1*e1*e3 - S2*e1*e2 + S2*e3 + S3*e1**2 - S3*e2
S_6 = S1*e1**2*e3 - S1*e2*e3 - S2*e1**2*e2 + S2*e1*e3 + S2*e2**2 + S3*e1**3 - 2*S3*e1*e2 + S3*e3

.........

Verifica usando le radici reali

Dopo aver calcolato le somme di potenze tramite le identità di Newton, è naturale chiedersi:

“Il risultato è davvero corretto?”

Per verificarlo possiamo usare un approccio completamente diverso:

1. trovare numericamente le radici del polinomio
2. elevare direttamente le radici alla potenza desiderata
3. sommare i risultati
4. confrontare con la formula ricorsiva

Questa è una verifica numerica indipendente.

Troviamo le radici

roots = sp.nroots(P)

print(roots)

[-0.801937735804838, 0.554958132087371, 2.24697960371747]

Cosa fa nroots()?

La funzione:

sp.nroots()

significa:

numeric roots

cioè:

“trova approssimazioni numeriche delle radici”.

Differenza con solve()

solve()

cerca formule esatte:

[math]\displaystyle \begin{aligned}
\sqrt{2}, \quad \frac{1+\sqrt{5}}{2}, \dots
\end{aligned}[/math]

nroots()

calcola approssimazioni decimali:

1.7320508,  0.6180339,…

Perché usare nroots()?

Perché molti polinomi:

• non hanno formule semplici
• hanno radici irrazionali
• hanno radici complesse

nroots() usa metodi numerici iterativi molto efficienti.

Verifica numerica

S5_numeric = sum(r**5 for r in roots)

print(sp.N(S5_numeric))

57.0000000000000

57.0000000000000