Minimi quadrati.

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica Pura e Applicata

Esempio 1

21 novembre 2011 Nel 1601 Keplero formulò la terza legge del moto planetario

T = Cx^{3/2}

dove x è la distanza dal sole in milione di km., T è il periodo dell’orbita misurato in giorni e C è una costante. Le coppie (x, T ) osservate per gli 8 planeti e Plutone sono

Mercurio 57.59 87.99
Venere 108.11 224.70
Terra 149.57 365.26
Marte 227.84 686.98
Giove 778.14 4332.4
Saturno 1427.0 10759
Urano 2870.3 30684
Nettuno 4499.9 60188
Plutone 5909.0 90719

Determinare C cosicchè T = Cx^{3/2}.

Alvise Sommariva  Minimi quadrati.  2/ 33

Esempio 1

La costante C non è immediatamente calcolabile. Infatti:

>> x = [ 5 7 . 5 9 ; 1 0 8 . 1 1 ; 1 4 9 . 5 7 ; 2 2 7 . 8 4 ; 7 7 8 . 1 4 ; …

1 4 2 7 . 0 ; 2 8 7 0 . 3 ; 4 4 9 9 . 9 ; 5 9 0 9 . 0 ] ;

>> T = [ 8 7 . 9 9 ; 2 2 4 . 7 0 ; 3 6 5 . 2 6 ; 6 8 6 . 9 8 ; 4 3 3 2 . 4 ; …
10759 ;30684 ;60188 ;90719 ] ;
>> format long e
>> C=T . / ( x . ˆ ( 3 / 2 ) )
C = 2.013319499958850 e−01
1.998960771429300 e−01
1.996803707178690 e−01
1.997556793162881 e−01
1.995914513555481 e−01
1.995886992640107 e−01
1.995358969439935 e−01
1.993910306711459 e−01
1.997226267954996 e−01

>>

Qual’e’ il valore attribuire a C?

Alvise Sommariva  Minimi quadrati.  3/ 33

Siano m la magnitudine apparente di una cometa, Δ e r rispettivamente le distanze in unità astronomiche della cometa dalla Terra e dal Sole, g la magnitudine assoluta e ξ un coefficiente. Noto che

m = g + 5 log_{10} (Δ) + ξ log_{10} (r)

e determinati alcuni valori di m, Δ, r, determinare g e ξ. Quindi si desiderano calcolare i migliori A ≈ g, B ≈ ξ tali che

A + B log_{10} (r) ≈ m − 5 log_{10} (Δ).

Problema: cosa vuol dire migliori?

Esempio. La tabella seguente riporta le stime di magnitudine visuale m della cometa periodica Wild 2, rilevate da J. Bortle, nochè r e Δ. I valori corrispondenti sono stati calcolati dagli elementi orbitali (IAUC 3177). Calcolare g e ξ.

1978 FEB. FEB. FEB. FEB. FEB. MAR. MAR. MAR.
UT 4.01 5.00 9.02 10.02 25.03 7.07 14.03 30.05
m 11.4 11.5 11.5 11.3 11.5 11.3 11.5 11.0
r  1.987 1.981 1.958 1.952 1.865 1.809 1.772 1.693
Δ 1.249 1.252 1.266 1.270 1.335 1.382 1.415 1.487
1978 APR. APR. APR. MAG. MAG. MAG. GIU. GIU.
UT 3.05 10.06 26.07 1.08 8.07 26.09 2.10 6.09
m 11.1 10.9 10.7 10.6 10.7 10.8 10.5 10.4
r  1.674 1.643 1.582 1.566 1.545 1.507 1.498 1.495
Δ 1.504 1.532 1.592 1.610 1.634 1.696 1.721 1.736

Plot risultati

Figura: Grafico che illustra l’approssimazione ai minimi quadrati A + B log_{10} (r) ≈ m − 5 log_{10} (Δ). Ascissa: r, ordinata: curva e m − 5 log_{10} (Δ).

Approssimazione ai minimi quadrati

Questi due problemi rientrano nella famiglia più ampia che segue.

▶ Sia fissata una certa funzione f : Ω ⊂ ℝ → ℝ di cui è noto il valore nei punti {x_{k}}_{k=1,...,n} ⊂ Ω. Hp. f ∈ C(Ω), Ω = (a, b) anche non limitato (per semplicità , si può generalizzare).

▶ Siano date m funzioni lin. indip. ϕ_{1}, …, ϕ_{m} : Ω ⊂ ℝ → ℝ. Si cerchino a = (a_{j})_{j=1,...,m} per cui la funzione

ψ^{∗}(x) = ∑_{j=1}^{m}a_{j}ϕ_{j}(x)

minimizza tra tutte le funzioni del tipo ψ(x) = ∑_{j=1}^{m} c_{j}ϕ_{j}(x),

∥f −ψ∥_{2,d} := \sqrt{∑_{k=1}^{n}|f (x_{k}) − ψ(x_{k})|^{2}} = \sqrt{∑_{k=1}^{n}|f (x_{k}) − ∑_{j=1}^{m}c_{j}ϕ_{j}(x_{k})|^{2}}

La funzione ψ^{∗} si dice approssimante ai minimi quadrati (discreti)

di f nello spazio S = span({ϕ_{k}}_{k=1,...,m}).

Alvise Sommariva  Minimi quadrati.  7/ 33

Approssimazione ai minimi quadrati e sistemi lineari

Posto ∥x∥_{2} := \sqrt{∑_{k=1}^{n} |x_{k}|^{2}} e definiti ▶ V la matrice n × m le cui componenti sono V_{i,j} = ϕ_{j}(x_{i}), ▶ y il vettore n × 1 le cui componenti sono y_{k} = f (x_{k}), ▶ a il vettore m × 1 le cui componenti sono a_{k},

il problema di approssimazione ai minimi quadrati consiste nel risolvere il sistema sovradeterminato V a = y cosicchè sia minima ∥V a − y∥_{2} in quanto

∥V a − y∥_{2} = \sqrt{∑_{k=1}^{n}|y_{k} − (V a)_{k}|^{2}} = \sqrt{∑_{k=1}^{n}|y_{k} − ∑_{j=1}^{m}V_{k,j}a_{j}|^{2}}

= \sqrt{∑_{k=1}^{n}|f (x_{k}) − ∑_{j=1}^{m}a_{j}ϕ_{j}(x_{k})|^{2}} = ∥f − ∑_{j=1}^{m}a_{j}ϕ_{j}∥_{2,d}

NB: Si osservi che V è rettangolare e il problema V a = y potrebbe non avere sol. classica a.

Alvise Sommariva  Minimi quadrati.  8/ 33

Funzione peso e prodotto scalare

Al momento non risulta chiaro come risolvere questo sistema sovradeterminato V a = y cosicchè sia minima ∥V a − y∥_{2}.

Sia w : Ω = (a, b) → ℝ una funzione peso (continua, strettamente positiva in (a, b), singolarità integrabili agli estremi sono ammesse).

▶ Jacobi: w(x) = (1 − x)^{α}(1 + x)^{β}, x ∈ (−1, 1), α, β > −1. Casi speciali: α = β = 0 (Legendre), α = β = −1/2 (Chebyshev), α = β (Gegenbauer).

▶ Gauss: w(x) = exp(−x^{2}), x ∈ ℝ.

Due classici prodotti scalari di due funzioni w-integrabili sono ▶ (f , g)_{c} := ∫_{Ω} f (x)g(x)w(x)dx, caso continuo; ▶ (f , g)_{d} := ∑_{i=1}^{n} f (x_{i})g(x_{i})w_{i}, w_{i} > 0, caso discreto.

I coeffs. w_{i} del caso discreto dipendono dalla funz. peso w. Proprietá: commutativa, linearitá, (f , f ) ≥ 0 per ogni f .

Alvise Sommariva  Minimi quadrati.  9/ 33

Norme, prod. scalari e ortogonalità

Definiamo ∥f ∥ = \sqrt{(f , f )}. Osserviamo che:

▶ ∥x∥_{2} := \sqrt{∑_{k=1}^{n} |x_{k}|^{2}} = \sqrt{(x, x)_{ℝ^{n}}} dove (·, ·)_{ℝ^{n}} è l’usuale

prodotto scalare di ℝ^{n}. ▶ ∥f − ψ∥_{2,d} := \sqrt{∑_{k=1}^{n} |f (x_{k}) − ψ(x_{k})|^{2}} = \sqrt{(f − ψ, f − ψ)_{d}}

Alcune definizioni ulteriori.

▶ Due funzioni f , g si dicono ortogonali se (f , g) = 0. ▶ Una sequenza finita o infinita di funzioni {ϕ_{k}}_{k=1,...} si dice

sistema ortogonale se (ϕ_{i}, ϕ_{j}) = 0 per i ≠ j, (ϕ_{i}, ϕ_{i}) > 0. ▶ Un sistema ortogonale {ϕ_{k}}_{k=1,...} per cui (ϕ_{i}, ϕ_{i}) = 1 si dice

sistema ortonormale. ▶ Una base che in particolare forma un sistema ortogonale

(ortonormale) si dice base ortogonale (base ortonormale).

Riassunto

▶ Dobbiamo cercare una funzione ψ^{∗} = ∑_{k=1}^{m} a_{k}ϕ_{k} che meglio approssimi f nei punti 2 a 2 distinti {x_{k}}, nel senso di minimizzare

∥f − ψ∥_{2,d} := \sqrt{∑_{k=1}^{n}|f (x_{k}) − ψ(x_{k})|^{2}} = \sqrt{(f − ψ, f − ψ)_{d}},

con (f , g)_{d} := ∑_{i=1}^{n} f (x_{i})g(x_{i}), quindi per w_{i} = 1. ▶ Il problema si converte in una questione di algebra lineare:

risolvere il sistema sovradeterminato V a = y cosicchè sia minima

∥V a − y∥_{2} = \sqrt{(V a − y, V a − y)_{ℝ^{n}}} .

Alvise Sommariva  Minimi quadrati.  11/ 33

Spazi euclidei

Uno spazio vettoriale E si dice euclideo, se è dotato di prodotto scalare. Ricordiamo che un prod. scalare è una funzione reale definita sulle coppie x, y ∈ E con le seguenti proprietà

1. (x, x) ≥ 0 per ogni x ∈ E; inoltre (x, x) = 0 se e solo se x = 0;

2. (x, y) = (y, x) per ogni x, y ∈ E;

3. (λx, y) = λ(x, y) per ogni x, y ∈ E e λ ∈ ℝ;

4. (x, y + z) = (x, y) + (x, z) per ogni x, y, z ∈ E.

Spazi euclidei

Vediamo alcuni notevoli esempi di spazio euclideo:

▶ Lo spazio ℝ^{m} dotato del prod. scal. (u, v)_{ℝ^{m}} := ∑_{k=1}^{m} u_{k} · v_{k}

è uno spazio euclideo. ▶ Sia Ω intv. compatto. L’insieme C(Ω) dotato del prodotto

scalare discreto (f , g)_{d} := ∑_{k=1}^{m} f (x_{k})g(x_{k}), con {x_{k}} ⊂ Ω, o del prod. scalare (f , g)_{c} := ∫_{Ω} f (x)g(x)dx è uno spazio euclideo.

▶ Sia Ω ⊆ ℝ un intervallo (anche illimitato), w : Ω → ℝ funz. peso. Sia L_{w}^{2}(Ω) lo spazio delle funzioni f : Ω ⊆ ℝ tali che ∫_{Ω} f^{2}(x)w(x)dx < +∞. L’insieme L_{w}^{2}(Ω) dotato del prodotto scalare continuo (f , g)_{c} := ∫_{Ω} f (x)g(x)w(x)dx è uno spazio euclideo.

Alcuni teoremi

Teorema. [Pitagora] Se (f , g) = 0 allora

∥f + g∥^{2} = ∥f ∥^{2} + ∥g∥^{2}.

Teorema. Se ϕ_{1}, …, ϕ_{n} è un sistema ortogonale, allora i ϕ_{j} sono linearmente indipendenti. Teorema. Se ϕ_{1}, …, ϕ_{n} è un sistema ortogonale, allora

∥ ∑_{j=1}^{m}c_{j}ϕ_{j}∥^{2} = (∑_{j=1}^{m}c_{j}ϕ_{j}, ∑_{k=1}^{m}c_{k}ϕ_{k}) = ∑_{j=1}^{m}c_{j}^{2}∥ϕ_{j}∥^{2}

per ogni {c_{j}}_{j=1,...,m}.

Alvise Sommariva  Minimi quadrati.  14/ 33

Teorema miglior approssimazione minimi quadrati

Teorema. Siano {ϕ_{j}}_{1,...,m} ∈ E funzioni linearmente indipendenti di uno spazio euclideo E. Sia f ∈ E. La sol. ψ^{∗} del problema

∥f − ψ^{∗}∥_{2} = min_{ψ∈span{ϕ_{j}}_{j=1,...,m}} ∥f − ψ∥_{2}

è

ψ^{∗} = ∑_{1,...,m}c_{j}^{∗}ϕ_{j}

dove i coefficienti c_{j}^{∗} verificano le cosidette equazioni normali

∑_{k=1}^{m}(ϕ_{j}, ϕ_{k})c_{k}^{∗} = (ϕ_{j}, f ), j = 1, …, m.  (1)

La soluzione è caratterizzata dalla proprietà di ortogonalitá cioè che (ψ^{∗} − f , ϕ_{k}) = 0, con k = 1, …, m, o equivalentemente

(ψ^{∗}, ϕ_{k}) = (f , ϕ_{k}), k = 1, …, m.  (2)

Teorema miglior approssimazione minimi quadrati: dim.

Supponiamo ψ^{∗} = ∑_{1,...,m} c_{j}^{∗}ϕ_{j} verifichi (ψ^{∗} − f , ϕ_{k}) = 0, con k = 1, …, m. Sia c = (c_{k})_{1,...,m}, una vettore di coefficienti e supponiamo che per almeno un indice j sia c_{j} ≠ c_{j}^{∗}, cioè c ≠ c^{∗} = (c_{k}^{∗})_{1,...,m}. Allora

∑_{j=1}^{m}c_{j}ϕ_{j} − f = (∑_{j=1}^{m}c_{j}ϕ_{j} − ψ^{∗}) + (ψ^{∗} − f )

= ∑_{j=1}^{m}(c_{j} − c_{j}^{∗})ϕ_{j} + (ψ^{∗} − f )  (3)

Se u = ψ^{∗} − f è ortogonale a tutti i ϕ_{j}, allora è ortogonale pure alla combinazione lineare di ϕ_{j} come ad esempio

v = ∑_{j=1}^{m}(c_{j} − c_{j}^{∗})ϕ_{j} = ∑_{j=1}^{m}c_{j}ϕ_{j} − ψ^{∗} ∈ span {ϕ_{k}}_{k=1,...,m}.

Dal teorema di Pitagora, poichè (u, v) = 0 implica ∥u + v∥^{2} = ∥u∥^{2} + ∥v∥^{2}, da ∥ ∑_{j=1}^{m}(c_{j} − c_{j}^{∗})ϕ_{j}∥ > 0 poichè c ≠ c^{∗} abbiamo

∥ ∑_{j=1}^{m}c_{j}ϕ_{j} − f ∥^{2} = ∥ (∑_{j=1}^{m}c_{j}ϕ_{j} − ψ^{∗}) + (ψ^{∗} − f )∥^{2}

= ∥ ∑_{j=1}^{m}c_{j}ϕ_{j} − ψ^{∗}∥^{2} + ∥ψ^{∗} − f ∥^{2}

= ∥ ∑_{j=1}^{m}(c_{j} − c_{j}^{∗})ϕ_{j}∥^{2} + ∥ψ^{∗} − f ∥^{2}
> ∥ψ^{∗} − f ∥^{2}  (4)

Di conseguenza se ψ^{∗} ∈ span{ϕ_{k}}_{k=1,...,m} e ψ^{∗} − f è ortogonale a tutti i ϕ_{k} allora ψ^{∗} è la miglior approssimazione di f in span{ϕ_{k}}_{k=1,...,m}. Rimane allora da mostrare che le condizioni di ortogonalità

(∑_{j=1}^{m}c_{j}^{∗}ϕ_{j} − f , ϕ_{k} ) = 0, k = 1, …, m

possano essere soddisfatte per un qualche c^{∗} = (c_{j})_{j=1,...,m}. Questo problema è equivalente alla soluzione del sistema di equazioni normali

∑_{k=1}^{m}(ϕ_{j}, ϕ_{k})c_{k}^{∗} = (ϕ_{j}, f ), j = 1, …, m  (5)

Se ϕ_{1}, …, ϕ_{m} sono m vettori linearmente indipendenti che formano un sistema ortogonale, si ha che

∑_{k=1}^{m}(ϕ_{j}, ϕ_{k})c_{k}^{∗} = (ϕ_{j}, ϕ_{j})c_{j}^{∗}

e quindi da (5) che

(ϕ_{k}, ϕ_{k})c_{k}^{∗} = (ϕ_{k}, f ).

Visto che (ϕ_{k}, ϕ_{k}) ≠ 0 (se così non fosse 0 = (ϕ_{k}, ϕ_{k}) = ∥ϕ_{k}∥^{2} avremmo ϕ_{k} = 0 e quindi {ϕ_{j}}_{j=1,...,m} non sarebbe un sistema di vettori linearmente indipendenti) si vede subito che (c_{k}^{∗})_{k} esistono unici e uguali a

c_{k}^{∗} = \frac{(ϕ_{k}, f )}{(ϕ_{k}, ϕ_{k})}.

Se invece ϕ_{1}, …, ϕ_{m} non formano un sistema ortogonale, il sistema di equazioni normali ha una e una sola soluzione se il sistema omogeneo di equazioni

∑_{k=1}^{m}(ϕ_{j}, ϕ_{k})c_{k}^{∗} = 0, j = 1, …, m  (6)

ha la sola soluzione nulla.

Teorema miglior approssimazione minimi quadrati: dim.

Se così non fosse, esisterebbe c^{∗} = (c_{j})_{j=1,...,m} per cui da (6)

∥ ∑_{j=1}^{m}c_{j}ϕ_{j}∥^{2} = (∑_{j=1}^{m}c_{j}ϕ_{j}, ∑_{k=1}^{m}c_{k}ϕ_{k} )
= ∑_{k=1}^{m}c_{k} ∑_{j=1}^{m}c_{j} (ϕ_{j}, ϕ_{k})
= ∑_{k=1}^{m}c_{k} · 0
= 0  (7)

e quindi essendo ∥ · ∥ una norma, necessariamente ∑_{j=1}^{m} c_{j}^{∗}ϕ_{j} = 0, il che contraddice il fatto che i ϕ_{k} erano linearmente indipendenti.

Alvise Sommariva  Minimi quadrati.  21/ 33

Corollario miglior approssimazione minimi quadrati

Teorema. Se {ϕ_{j}}_{j=1,...,m} è un sistema ortogonale allora i coefficienti c_{j}^{∗} (detti in questo caso di Fourier) sono calcolabili più semplicemente con la formula

c_{j}^{∗} = \frac{(f , ϕ_{j})}{(ϕ_{j}, ϕ_{j})}, j = 1, …, m.

Risoluzione minimi quadrati iniziale

▶ Si determina una base ortogonale {ϕ_{k}^{⊥}} dello spazio vett. S = span({ϕ_{k}}_{k=1,...,m}) (algoritmo di Gram-Schmidt generalizzato);

▶ Si calcola c_{j}^{∗} = \frac{(f ,ϕ_{j})}{(ϕ_{j},ϕ_{j})}, j = 1, …, m;

▶ Si ottiene ψ^{∗} = ∑_{j} c_{j}^{∗}ϕ_{j}^{⊥}.

Minimi quadrati in Matlab/Octave

Per risolvere il problema ai minimi quadrati in Matlab (sistema sovradeterminato), basta usare il comando \.

Risolviamo l’esempio 1, salvando nel file esempio1.m il seguente codice.

x = [ 5 7 . 5 9 ; 1 0 8 . 1 1 ; 1 4 9 . 5 7 ; 2 2 7 . 8 4 ; 7 7 8 . 1 4 ; …

1 4 2 7 . 0 ; 2 8 7 0 . 3 ; 4 4 9 9 . 9 ; 5 9 0 9 . 0 ] ;

T = [ 8 7 . 9 9 ; 2 2 4 . 7 0 ; 3 6 5 . 2 6 ; 6 8 6 . 9 8 ; 4 3 3 2 . 4 ; …
10759 ;30684 ;60188 ;90719 ] ;
%T=C x ˆ {3/2}.
V=x . ˆ ( 3 / 2 ) ; C=V\T ;
format long e ;

C

[ T C∗V ]

Risultati esempio 2

% d a t i =[UT m r d e l t a ]
dati=[ 4.01  11.4  1.987  1.249
5.00  11.5  1.981  1.252
9.02  11.5  1.958  1.266
10.02 11.3  1.952  1.270
25.03 11.5  1.865  1.335
7.07  11.3  1.809  1.382
14.03 11.5  1.772  1.415
30.05 11.0  1.693  1.487
3.05  11.1  1.674  1.504
10.06 10.9  1.643  1.532
26.07 10.7  1.582  1.592
1.08  10.6  1.566  1.610
8.07  10.7  1.545  1.634
26.09 10.8  1.507  1.696
2.10  10.5  1.498  1.721
6.09  10.4  1.495  1 . 7 3 6 ] ;
UT=dati ( : , 1 ) ; m=dati ( : , 2 ) ; r=dati ( : , 3 ) ; delta=dati ( : , 4 ) ;
V=[ones ( s i z e ( UT ) ) l o g 1 0 ( r ) ] ; y=m−5∗l o g 1 0 ( delta ) ; c=V\y ;
format long e ; c , [ m−5∗l o g 1 0 ( delta ) V∗c ]
Alvise Sommariva  Minimi quadrati.  25/ 33

Risultati esempio 2

>> esempio2

c = 7.065569523839867 e+00
1.346222475779441 e+01
ans = 1.091718780812932 e+01  1.107997623316711 e+01
1.101197835562795 e+01  1.106229506063970 e+01
1.098783147159332 e+01  1.099401750969824 e+01
1.078098139522022 e+01  1.097607405849955 e+01
1.087259367149703 e+01  1.070950885319811 e+01
1.059745978480910 e+01  1.053126537223116 e+01
1.074621780069845 e+01  1.041044393366901 e+01
1.013844515739023 e+01  1.014380088633594 e+01
1.021376081872188 e+01  1.007781573175590 e+01
9.973706173517074 e+00  9.968530869025818 e+00
9.690284682991749 e+00  9.747331919522992 e+00
9.565870619840751 e+00  9.687899947258250 e+00
9.633739739018015 e+00  9.608967235377786 e+00
9.652870760396526 e+00  9.463370158841187 e+00
9.321095648362199 e+00  9.428349002986664 e+00
9.202251395797635 e+00  9.416628511404101 e+00

>>

Alvise Sommariva  Minimi quadrati.  26/ 33

Minimi quadrati e smoothing

Digitiamo sulla shell di Matlab

>> x =0:0.01:2∗ p i ;
>> y=s i n (2∗x ) +(10ˆ (−1) ) ∗rand ( s i z e ( x ) ) ;
>> p l o t ( x , y , ’ r−’ ) ;

▶ Interpretazione: perturbazione della funzione sin(2x)

nell’intervallo [0, 2π]. ▶ Necessitá: ricostruire sin(2x) (e non funz. perturbata). ▶ Nota: non ha senso utilizzare un interpolante polinomiale p di

grado n nè una spline interpolante visto che ricostruirebbero la funzione perturbata.

Minimi quadrati e polyfit

Scriviamo sulla shell di Matlab/Octave help polyfit. In una recente release di Matlab appare

 POLYFIT Fit polynomial to data .
     POLYFIT ( X , Y , N ) finds the coefficients of a
         polynomial P ( X )
     of degree N that fits the data , P ( X ( I ) )˜ =Y ( I ) , in a

          least

    −squares sense .

In altri termini polyfit calcola i coefficienti del polinomio p_{N} di grado N che meglio approssima (in norma 2 discreta) la funzione f avente nel vettore di nodi X i valori Y (cioè Y (i) := f (X (i))). Operativamente si cerca il polinomio p_{N} per cui risulta minima

∥f − p_{N}∥_{2,d} = \sqrt{∑_{i}|f (x_{i}) − p_{N}(x_{i})|^{2}}.

Alvise Sommariva  Minimi quadrati.  28/ 33

Applicazione polyfit

▶ Posto f (x) = sin (x) e f̃(x) = sin (x) + δ(x), confrontiamo graficamente per n = 2, …, 8, nei nodi x_{k} = kh, h = 2π/999, k = 0, …, 999, la funzione perturbata f (x) = sin (x) + δ(x) con la approssimante ai minimi quadrati p_{n}^{∗}.

▶ Valutiamo per n = 2, …, 8 la quantità ∥f − p_{n}^{∗}∥_{2,d} e

∥f̃ − p_{n}^{∗}∥_{2,d}. ▶ Valutiamo per n = 2, …, 8 la quantità ∥f − p_{n}^{∗}∥_{∞,d} e

∥f̃ − p_{n}^{∗}∥_{∞,d}, dove ricordiamo

∥g∥_{∞,d} = max_{k=0,...,999}|g(x_{k})|.

Applicazione polyfit

Salviamo in esempio3.m:

x=l i n s p a c e (0 ,2∗ pi , 1 0 0 0 ) ;
y=s i n (2∗x ) ; yy=y+(10ˆ (−2) ) ∗rand ( s i z e ( x ) ) ;
f o r n=2:8
coeff=p o l y f i t ( x , yy , n ) ; % COEFFS . BEST APPROX (B . A . )
z=p o l y v a l ( coeff , x ) ;  % VALORE B . A . NEI NODI ”x ” .
p l o t ( x , yy , ’ r−’ , x , z , ’ k−’ ) ;
err2=norm ( z−y , 2 ) ; err2p=norm ( z−yy , 2 ) ; % ERRS .

errinf=norm ( z−y , inf ) ; errinfp=norm ( z−yy , inf ) ; %

ERRS .
f p r i n t f ( ’ \n\ t [DEG] : % 2 . 0 f ’ , n ) ;
f p r i n t f ( ’ [ 2 ] : % 2 . 2 e %2.2 e ’ , err2 , err2p ) ;
f p r i n t f ( ’ [ INF ] : % 2 . 2 e %2.2 e ’ , errinf , errinfp ) ;
pause ( 2 ) ;

end

f p r i n t f ( ’ \n \n ’ ) ;

Risultati

>> % [ ERR . ] [ 2 ] : SIN.−SIN PERT . [ INF ] : SIN.−SIN PERT .

>> esempio3

[ DEG ] : 2 [ 2 ] : 2 . 0 6 e+01 2.06 e+01 [ INF ] : 1 . 1 3 e+00 1.13 e+00
[ DEG ] : 3 [ 2 ] : 1 . 8 9 e+01 1.89 e+01 [ INF ] : 1 . 1 6 e+00 1.16 e+00
[ DEG ] : 4 [ 2 ] : 1 . 8 9 e+01 1.89 e+01 [ INF ] : 1 . 1 6 e+00 1.16 e+00
[ DEG ] : 5 [ 2 ] : 6 . 8 6 e+00 6.86 e+00 [ INF ] : 6 . 6 7 e−01 6.64 e−01
[ DEG ] : 6 [ 2 ] : 6 . 8 6 e+00 6.86 e+00 [ INF ] : 6 . 6 8 e−01 6.64 e−01
[ DEG ] : 7 [ 2 ] : 1 . 2 3 e+00 1.22 e+00 [ INF ] : 1 . 4 8 e−01 1.44 e−01
[ DEG ] : 8 [ 2 ] : 1 . 2 3 e+00 1.22 e+00 [ INF ] : 1 . 4 7 e−01 1.44 e−01

>>

Plot risultati

Figura: Grafico che illustra l’approssimazione ai minimi quadrati di grado 5 su una perturbazione della funzione sin (2x) (campionamento in nodi equispaziati)).

Nota: polyfit e interpolazione

Supponiamo fissati i punti {x_{k}}_{k=1,...,m} (a due a due distinti) e sia p_{m−1} il polinomio che interpola le coppie (x_{k},f (x_{k})) per k = 1, …, m Evidentemente da f (x_{i}) = p_{m−1}(x_{i}) per i = 1, …, m abbiamo

∥f − p_{m−1}∥_{2,d} = \sqrt{∑_{i=1}^{m}|f (x_{i}) − p_{m−1}(x_{i})|^{2}} = 0,

e quindi il polinomio interpolatore risulta la approssimante ai minimi quadrati di f (relativa alla norma 2 discreta basata sui punti {x_{k}}_{k=1,...,m}).

Di conseguenza, il comando coeffs=polyfit(x,y,n-1) darà i coefficienti del polinomio interpolatore qualora i vettori x, y abbiano dimensione n.