Appendice A — Complementi sui polinomi ortogonali
Questa appendice raccoglie quattro risultati sui polinomi ortogonali che non servono per i minimi quadrati di Capitolo 17 né per le quadrature gaussiane di Capitolo 18, e che sono stati tolti dal Capitolo 16 per non appesantirlo. Restano qui perché sono strumenti di calcolo standard: si incontrano in ogni trattazione dell’argomento e servono appena si va un passo oltre il programma.
A.1 Simmetria: cosa succede se il peso è pari
Supponiamo che l’intervallo sia simmetrico rispetto all’origine, (a,b)=(-c,c), e che il peso sia una funzione pari,
w(-x)=w(x)\qquad\text{per ogni }x .
È il caso dei polinomi di Legendre (w\equiv1 su [-1,1]), di Chebyshev (w=(1-x^2)^{-1/2}) e in generale di Gegenbauer. Sotto questa ipotesi i polinomi ortogonali sono alternativamente pari e dispari:
p_n(-x)=(-1)^n\,p_n(x). \tag{A.1}
Vediamo perché, usando i polinomi monici \pi_n e l’unicità garantita da Teorema 14.1. Definiamo il polinomio «riflesso»
q_n(x)=(-1)^n\,\pi_n(-x)
e verifichiamo che possiede le due proprietà che caratterizzano \pi_n.
È monico di grado n. Sostituire x\mapsto-x in \pi_n moltiplica il termine di testa x^n per (-1)^n; il fattore (-1)^n davanti lo riporta a +1.
È ortogonale a tutti i polinomi di grado minore. Sia r con \deg r<n. Nel calcolo
\langle q_n,r\rangle_w=(-1)^n\int_{-c}^{c}\pi_n(-x)\,r(x)\,w(x)\,\mathrm dx
operiamo il cambio di variabile t=-x. L’intervallo si ribalta in sé (è simmetrico), \mathrm dx diventa -\mathrm dt con scambio degli estremi, e w(-t)=w(t) perché il peso è pari:
\langle q_n,r\rangle_w=(-1)^n\int_{-c}^{c}\pi_n(t)\,r(-t)\,w(t)\,\mathrm dt =(-1)^n\big\langle \pi_n,\,r(-\,\cdot\,)\big\rangle_w=0,
perché anche r(-\,\cdot\,) ha grado minore di n e \pi_n è ortogonale a tutti questi polinomi. Dunque q_n è monico di grado n e ortogonale a \mathcal P_{n-1}: per l’unicità q_n=\pi_n, che è esattamente A.1.
Da A.1 si leggono tre conseguenze immediate.
- Gli zeri sono simmetrici. Se p_n(x_0)=0 allora anche p_n(-x_0)=\pm p_n(x_0)=0: gli zeri si presentano in coppie \pm x_0.
- Per n dispari, lo zero è x=0. Ponendo x=0 in A.1 con n dispari si ottiene p_n(0)=-p_n(0), cioè p_n(0)=0. L’origine è dunque sempre un nodo delle formule gaussiane di ordine dispari.
- Nella ricorrenza b_k=0. Infatti b_k=\langle x\phi_k,\phi_k\rangle_w è l’integrale di x\,\phi_k(x)^2w(x), prodotto di una funzione dispari (x) per una pari (\phi_k^2w): l’integrale su un intervallo simmetrico di una funzione dispari è nullo. Si ritrova così la 16.7: la matrice di Jacobi dei polinomi di Legendre ha diagonale nulla.
Per i Jacobi con \alpha\ne\beta il peso non è pari, e infatti la simmetria non vale nella forma A.1. La sostituzione x\mapsto-x trasforma però il peso in
(1+x)^\alpha(1-x)^\beta=w_{\beta,\alpha}(x),
cioè nel peso con i due parametri scambiati. Ripetendo l’argomento di sopra si ottiene
P_n^{(\alpha,\beta)}(-x)=(-1)^n\,P_n^{(\beta,\alpha)}(x). \tag{A.2}
Per \alpha=\beta i due membri riguardano la stessa famiglia e si ricade in A.1.
A.2 Derivate dei polinomi ortogonali
Serve spesso la derivata di un polinomio ortogonale: per esempio nella formula dei pesi di Gauss–Legendre \omega_i=2/\big[(1-x_i^2)P_n'(x_i)^2\big] di Capitolo 18, o per raffinare gli zeri con il metodo di Newton. Derivare la forma monomiale sarebbe disastroso; conviene una formula che riusi i valori già calcolati con la ricorrenza.
Per i polinomi di Legendre vale
(1-x^2)\,P_n'(x)=n\big[P_{n-1}(x)-x\,P_n(x)\big], \qquad -1<x<1 , \tag{A.3}
che si ricava dalla formula di Rodrigues 14.6 con alcune integrazioni per parti (il conto è in [1]). Dal punto di vista pratico A.3 dice che, una volta ottenuti P_{n-1}(x) e P_n(x) con la ricorrenza a tre termini, la derivata costa due moltiplicazioni e una divisione, senza alcun lavoro aggiuntivo.
Agli estremi x=\pm1 la formula va letta come un limite, perché il fattore 1-x^2 si annulla; il valore corretto è
P_n'(1)=\frac{n(n+1)}{2}, \qquad P_n'(-1)=(-1)^{n-1}\frac{n(n+1)}{2}.
Per i Jacobi vale invece una formula di natura diversa, che non lega la derivata alla stessa famiglia ma a una famiglia vicina:
\frac{\mathrm d}{\mathrm dx}P_n^{(\alpha,\beta)}(x) =\frac{n+\alpha+\beta+1}{2}\,P_{n-1}^{(\alpha+1,\beta+1)}(x). \tag{A.4}
Vale la pena fermarsi un attimo su cosa dice. Derivando si abbassa il grado di uno, come sempre; ma il risultato non è un polinomio qualunque: è di nuovo un polinomio di Jacobi, con entrambi i parametri aumentati di uno, cioè ortogonale rispetto al peso (1-x)^{\alpha+1}(1+x)^{\beta+1}, che si annulla agli estremi più rapidamente. In altre parole la derivazione manda famiglie ortogonali in famiglie ortogonali: una proprietà tutt’altro che scontata, che è alla base delle formule di quadratura di Gauss–Lobatto e dei metodi spettrali.
A.3 La formula di Christoffel–Darboux
Nel Sezione 16.3 si incontra ripetutamente la somma
\sum_{k=0}^{n-1}\frac{p_k(x)\,p_k(y)}{h_k}, \qquad h_k=\langle p_k,p_k\rangle_w , \tag{A.5}
dove p_k è la famiglia monica. Calcolarla come sta scritta richiede di valutare n polinomi. La formula di Christoffel–Darboux12 riduce tutto a due valutazioni:
\sum_{k=0}^{n-1}\frac{p_k(x)p_k(y)}{h_k}= \frac{p_n(x)\,p_{n-1}(y)-p_{n-1}(x)\,p_n(y)}{h_{n-1}\,(x-y)}, \qquad x\ne y . \tag{A.6}
La dimostrazione è un conto diretto, che conviene seguire passo passo.
Passo 1: la ricorrenza monica e un’identità sui h_k. Nella normalizzazione monica la ricorrenza si scrive
p_{k+1}(t)=(t-b_k)\,p_k(t)-a_k^2\,p_{k-1}(t), \qquad p_{-1}=0,\ p_0=1 . \tag{A.7}
Ci servirà il valore di a_k^2. Moltiplichiamo scalarmente A.7 per p_{k-1}: a sinistra si ottiene 0 perché p_{k+1}\perp p_{k-1}, e a destra sopravvivono due termini,
0=\langle t\,p_k,p_{k-1}\rangle_w-a_k^2\,h_{k-1} .
Ora \langle t\,p_k,p_{k-1}\rangle_w=\langle p_k,\,t\,p_{k-1}\rangle_w e t\,p_{k-1} è monico di grado k, dunque uguale a p_k più un polinomio di grado minore, che è ortogonale a p_k. Resta \langle t\,p_k,p_{k-1}\rangle_w=h_k e quindi
a_k^2=\frac{h_k}{h_{k-1}} . \tag{A.8}
Passo 2: scrivere la ricorrenza in due punti e sottrarre. Scriviamo A.7 una volta in t=x e una volta in t=y:
\begin{aligned} p_{k+1}(x)&=(x-b_k)p_k(x)-a_k^2p_{k-1}(x),\\[2pt] p_{k+1}(y)&=(y-b_k)p_k(y)-a_k^2p_{k-1}(y). \end{aligned}
Moltiplichiamo la prima per p_k(y), la seconda per p_k(x) e sottraiamo. I termini con b_k compaiono come -b_kp_k(x)p_k(y) in entrambe le righe e si cancellano:
\underbrace{p_{k+1}(x)p_k(y)-p_k(x)p_{k+1}(y)}_{=:\,F_k} =(x-y)\,p_k(x)p_k(y) +a_k^2\underbrace{\big[p_k(x)p_{k-1}(y)-p_{k-1}(x)p_k(y)\big]}_{=\,F_{k-1}} . \tag{A.9}
Conviene dare un nome alla combinazione che si ripete: posto
F_k=p_{k+1}(x)\,p_k(y)-p_k(x)\,p_{k+1}(y),
la A.9 diventa semplicemente
F_k=(x-y)\,p_k(x)p_k(y)+a_k^2\,F_{k-1} .
Passo 3: dividere per h_k e far telescopare. Dividiamo per h_k e usiamo A.8, cioè a_k^2/h_k=1/h_{k-1}:
\frac{F_k}{h_k}-\frac{F_{k-1}}{h_{k-1}} =(x-y)\,\frac{p_k(x)p_k(y)}{h_k} . \tag{A.10}
A sinistra compare la differenza di due termini consecutivi della stessa successione k\mapsto F_k/h_k: sommando su k=0,1,\ldots,n-1 tutti i termini intermedi si elidono e resta solo la differenza fra l’ultimo e il primo,
\frac{F_{n-1}}{h_{n-1}}-\frac{F_{-1}}{h_{-1}} =(x-y)\sum_{k=0}^{n-1}\frac{p_k(x)p_k(y)}{h_k} .
Il termine F_{-1}=p_0(x)p_{-1}(y)-p_{-1}(x)p_0(y) è nullo perché p_{-1}=0. Esplicitando F_{n-1} e dividendo per x-y\ne0 si ottiene A.6.
Il caso y=x. Per y\to x numeratore e denominatore di A.6 tendono entrambi a zero. Applicando la regola di de l’Hôpital rispetto a y (o, equivalentemente, riconoscendo nel rapporto un rapporto incrementale) si ottiene la forma confluente
\sum_{k=0}^{n-1}\frac{p_k(x)^2}{h_k} =\frac{p_n'(x)\,p_{n-1}(x)-p_{n-1}'(x)\,p_n(x)}{h_{n-1}} , \tag{A.11}
Il membro sinistro è una somma di quadrati, quindi strettamente positivo: è la quantità che compare al denominatore della funzione di Christoffel 16.17 e, attraverso di essa, nei pesi delle formule gaussiane.
A.4 Completezza e proiezione ortogonale
Per un peso di Jacobi i polinomi sono densi in L^2_w(-1,1). La proiezione ortogonale di f\in L^2_w sulla classe \mathcal P_n è \Pi_nf=\sum_{k=0}^n c_k\phi_k,\qquad c_k=\langle f,\phi_k\rangle_w. Poiché f-\Pi_nf è ortogonale a \mathcal P_n, per qualunque q\in\mathcal P_n si ha la decomposizione di Pitagora \|f-q\|_w^2=\|f-\Pi_nf\|_w^2+\|\Pi_nf-q\|_w^2. Questa prova mostra che \Pi_nf è l’unica migliore approssimazione in norma pesata. La densità implica \|f-\Pi_nf\|_w\to0.
In sintesi
- Se l’intervallo è simmetrico e il peso è pari, i polinomi sono alternativamente pari e dispari: gli zeri vengono in coppie \pm x_0, l’origine è uno zero nei gradi dispari e la matrice di Jacobi ha diagonale nulla.
- Le derivate si ottengono senza derivare la forma monomiale: per i Legendre con A.3, per i Jacobi con A.4, che manda la famiglia (\alpha,\beta) nella famiglia (\alpha+1,\beta+1).
- La formula di Christoffel–Darboux riduce una somma di n termini a due sole valutazioni; la sua forma confluente dà il reciproco della funzione di Christoffel.
- La proiezione ortogonale è la migliore approssimazione in norma pesata, e la densità dei polinomi garantisce che l’errore tenda a zero.