ARGOMENTI DEL CORSO CALCOLO NUMERIC0vvw.web.cs.unibo.it/wiki/images/0/09/CC-2009_2010-02.pdf ·...

ARGOMENTI DEL CORSO

CALCOLO NUMERIC0

A.A. 2009/10

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Interpolazionecon Funzioni Polinomiali

Giulio Casciola

(ottobre 2003, rivista e corretta ottobre 2009)

Indice

1 Funzioni Polinomiali 11.1 Valutazione Numerica . . . . . . . . . . . . . . . . . . . . 3

1.1.1 Lo schema di Ruffini-Horner . . . . . . . . . . . . . 41.1.2 Valutazione della Derivata . . . . . . . . . . . . . . 5

1.2 Valutazione ed Errore Inerente . . . . . . . . . . . . . . . . 71.3 Polinomi di Bernstein . . . . . . . . . . . . . . . . . . . . . 10

1.3.1 Errore algoritmico . . . . . . . . . . . . . . . . . . 131.3.2 Proprieta . . . . . . . . . . . . . . . . . . . . . . . 151.3.3 Formula ricorrente . . . . . . . . . . . . . . . . . . 161.3.4 Valutazione numerica . . . . . . . . . . . . . . . . . 171.3.5 Suddivisione . . . . . . . . . . . . . . . . . . . . . . 191.3.6 Derivata . . . . . . . . . . . . . . . . . . . . . . . . 211.3.7 Antiderivata e integrazione . . . . . . . . . . . . . . 23

1.4 Un’applicazione: le curve di Bezier . . . . . . . . . . . . . 24

2 Interpolazione 272.1 Forma Monomiale . . . . . . . . . . . . . . . . . . . . . . . 282.2 Forma di Bernstein . . . . . . . . . . . . . . . . . . . . . . 292.3 Forma di Newton . . . . . . . . . . . . . . . . . . . . . . . 292.4 Forma di Lagrange . . . . . . . . . . . . . . . . . . . . . . 312.5 Errore nell’interpolazione polinomiale . . . . . . . . . . . . 322.6 Interpolazione polinomiale a tratti . . . . . . . . . . . . . . 35

2.6.1 Interpolazione locale . . . . . . . . . . . . . . . . . 362.6.2 Interpolazione global . . . . . . . . . . . . . . . . . 38

2.7 Un’applicazione: curve di interpolazione . . . . . . . . . . 44

3 Approssimazione ai minimi quadrati 473.1 Forma Monomiale . . . . . . . . . . . . . . . . . . . . . . . 493.2 Forma di Bernstein . . . . . . . . . . . . . . . . . . . . . . 493.3 Forma ortogonale . . . . . . . . . . . . . . . . . . . . . . . 50

3.3.1 Generazione polinomi ortogonali . . . . . . . . . . . 50

3

4 INDICE

3.4 Esempi numerici . . . . . . . . . . . . . . . . . . . . . . . 513.5 Retta di regressione lineare . . . . . . . . . . . . . . . . . . 533.6 Un’applicazione: curve di approssimazione . . . . . . . . . 54

4 Approssimazione di forma 554.1 Approssimazione uniforme . . . . . . . . . . . . . . . . . . 554.2 Approssimazione VD . . . . . . . . . . . . . . . . . . . . . 564.3 Significato geometrico dei coefficienti . . . . . . . . . . . . 58

Bibliografia 63

Capitolo 1

Funzioni Polinomiali nelCalcolo Numerico

Il problema dell’approssimazione di funzioni, cioe la determinazione diuna funzione piu semplice per rappresentarne una piu complessa, sele-zionandola da una prefissata classe di funzioni a dimensione finita, e unproblema di fondamentale importanza nella matematica applicata. Permeglio comprenderne l’importanza esaminiamo due situazioni classiche incui si presenta tale problema.

• La funzione da approssimare e nota in forma analitica, ma e richie-sta una sua approssimazione con una funzione piu semplice per lasua valutazione numerica, tipicamente con un calcolatore, oppureper renderne possibili operazioni quali ad esempio la derivazione ol’integrazione;

• La funzione f(x) da approssimare non e nota, ma di essa si cono-scono alcuni valori yi in corrispondenza di punti xi, i = 0, . . . , n e sivogliono avere indicazioni sul comportamento della funzione in altripunti. In questo caso la funzione e nota in senso discreto, per esem-pio un insieme di dati sperimentali, e approssimare la f(x) significadarne un modello rappresentativo.

In tal senso quello che vedremo in questa dispensa e un’introduzione allamodellazione. La Fig. 1.1 illustra un problema di interpolazione di puntinello spazio 3D con una funzione vettoriale (superficie in forma parame-trica) che solitamente e il primo passo in un processo di ricostruzione emodellazione di forma.

Il problema puo essere affrontato con differenti tecniche, che per essereadeguate devono tener conto della specifica situazione. In questo senso ri-sulta importante la scelta della distanza con la quale si vuole approssimare

1

2 CAPITOLO 1. FUNZIONI POLINOMIALI

Figura 1.1: Interpolazione di punti 3D con una funzione vettoriale

e che misura l’errore, mentre in certe situazioni, piuttosto che ridurre la di-stanza, risulta importante conservare certe proprieta di forma dei dati. Inquesta parte esamineremo la tecnica dell’interpolazione (o collocazione),sulla quale si basano la maggior parte dei metodi numerici per il calcolodegli integrali definiti e per l’approssimazione delle equazioni differenzia-li, la tecnica dell’approssimazione nel senso dei minimi quadrati e, menofrequente nei testi numerici, l’approssimazione di forma. Per tutte questetecniche ci limiteremo al caso polinomiale, ossia andremo a selezionarel’approssimazione desiderata dalla classe delle funzioni polinomiali a coef-ficienti reali e a variabile reale, nel semplice caso scalare. Questo capitolo,prima di introdurre le tecniche citate, vuole richiamare le nozioni piu im-portanti sulle funzioni polinomiali e discuterne il loro utilizzo dal puntodi vista numerico.

Definizione 1.1 (Polinomio) Una funzione p : IR→ IR definita da

p(x) = a0 + a1x + a2x2 + . . . + anxn =

n∑

i=0

aixi (1.1)

dove n e un intero non negativo e a0, a1, . . . , an sono numeri reali fissatie detto polinomio. In questa rappresentazione, se an 6= 0, si dice che p(x)ha grado n (ordine n + 1); se tutti i coefficienti ai sono nulli, allora p(x)e detto il polinomio nullo.

Con Pn si denota l’insieme di tutti i plinomi p con grado non superioread n, insieme al polinomio nullo. Si puo dimostrare che Pn e uno spaziovettoriale di dimensione n + 1.

La principale ragione dell’importanza delle funzioni polinomiali nelcalcolo numerico deriva dalle loro buone proprieta che qui richiamiamobrevemente:

1.1. VALUTAZIONE NUMERICA 3

• sono facilmente memorizzabili, e valutabili numericamente con uncalcolatore;

• sono funzioni regolari;

• la derivata e l’antiderivata di un polinomio sono ancora polinomi icui coefficienti sono determinati algoritmicamente (e quindi con uncalcolatore);

• il numero di zeri di un polinomio di grado n non puo essere superioread n;

• data una qualunque funzione continua su un intervallo [a, b], esisteun polinomio che la approssima uniformemente qualunque sia l’ǫprefissato

|f(x)− p(x)| < ǫ ∀x ∈ [a, b].

Teorema 1.1 (fondamentale dell’algebra) Sia p(x) un polinomio digrado n ≥ 1. L’equazione algebrica di grado n, p(x) = 0 ha almeno unaradice reale o complessa.

Corollario 1.1 Ogni equazione algebrica di grado n ha esattamente nradici reali o complesse, ciascuna con la sua molteplicita, cioe:

p(x) = an(x− α1)m1(x− α2)

m2 · · · (x− αk)mk (1.2)

dove αi, i = 1, . . . , k sono reali a due a due distinte e mi, i = 1, . . . , ksono le relative molteplicita, tali che m1 + m2 + · · ·+ mk = n.

Teorema 1.2 Siano a(x) e b(x) polinomi con b(x) 6= 0; allora esistono esono unici i polinomi q(x) ed r(x) per cui

a(x) = q(x)b(x) + r(x)

con r(x) ≡ 0 o r(x) con grado minore di quello di b(x).

1.1 Valutazione Numerica di un Polinomio

Assegnato un polinomio di grado n nella sua rappresentazione monomiale1.1, il metodo piu immediato per la sua valutazione in corrispondenza diun assegnato valore x puo essere descritto come segue:


s:=1

p:=a[0]

per k=1,..,n

s:=s*x_bar

p:=p+s*a[k]

p(x_bar):=p

Pertanto il calcolo di p(x) richiede 2n moltiplicazioni ed n addizio-ni/sottrazioni.Se si scrive il polinomio 1.1 nella seguente forma:

p(x) = a0 + x(a1 + x(a2 + · · ·+ x(an−1 + xan) · · ·)) (1.3)

si ricava il seguente metodo dovuto ad Horner:

p:=a[n]

per k=n-1,..,0

p:=a[k]+x_bar*p

p(x_bar):=p

L’algoritmo di Horner richiede n moltiplicazioni ed n addizioni/sottrazionie risulta quindi piu rapido del precedente oltre che numericamente piustabile.

1.1.1 Lo schema di Horner e la regola di Ruffini

In questa sezione si vuol ricavare il metodo di Horner da un contesto piugenerale che sara utile per ulteriori schemi di calcolo per polinomi. Richia-mando il Teorema 1.2 sulla divisione di due polinomi, nel caso particolarein cui il polinomio divisore sia il binomio (x− x) per un assegnato valorereale x, si ha che esistono unici i polinomi q(x) ed r(x) per cui

p(x) = q(x)(x− x) + r(x)

e poiche r(x) deve essere di grado inferiore a quello di (x − x) sara unacostante che indicheremo con r. Se si valuta p(x) in x si trova banalmenteche p(x) ≡ r. Quanto osservato viene formalizzato nel seguente teorema.

Teorema 1.3 Il resto della divisione del polinomio p(x) per (x − x) ep(x).

Dal teorema appena enunciato si deduce che un metodo per valutare unpolinomio p(x) in un punto x consiste nel determinare il resto della di-visione fra p(x) e il binomio (x − x). A tal fine e ben nota la regola diRuffini che viene applicata facendo uso del seguente schema di calcolo:

1.1. VALUTAZIONE NUMERICA 5

an an−1 ... ... a2 a1 a0

x xbn xbn−1 ... xb3 xb2 xb1

bn bn−1 ... ... b2 b1 r

In pratica, noti i coefficienti ai ed x, i bi ed r ≡ p(x) vengono ricavati nelseguente modo:

bn ← an

bn−1 ← an−1 + xbn

bn−2 ← an−2 + xbn−1...

......

b1 ← a1 + xb2

r ← a0 + xb1

dove

q(x) = b1 + b2x + b3x2 + . . . bnx

n−1 =n−1∑

i=0

bi+1xi.

Come si puo osservare, questo schema e esattamente quello di Horner visto;infatti se non interessa memorizzare i coefficienti del polinomio quozienteq(x) e li si sostituisce con una variabile semplice p i due schemi coincidono.

Esempio 1.1 Sia dato

p(x) = 1 + x− 2x2 + 3x4

e si voglia calcolare p(2). La regola di Ruffini e molto nota perche permettedi procedere anche manualmente in modo semplice;

3 0 −2 1 12 6 12 20 42

3 6 10 21 43 ≡ p(2)

1.1.2 Valutazione Numerica della Derivata

Si vuole ora procedere alla valutazione numerica della derivata di un po-linomio p(x) in un assegnato punto x, cioe si vuole valutare p′(x). Laprima idea potrebbe essere quella di calcolare numericamente i coefficientidel polinomio derivato e valutarlo in x con Horner; questo equivale a deter-minare l’espressione p′(x) = a1 +2a2x+3a3x

2 + . . . nanxn−1 e a procedere

nel valutare un polinomio di grado n − 1. Nel caso in cui oltre al valoredella derivata in x fosse necessario avere anche il valore del polinomio,allora, in modo piu economico, si puo procedere nel seguente modo: perquanto detto nella sezione precedente e

p(x) = q(x)(x− x) + p(x),


derivando tale espressione si ha

p′(x) = q′(x)(x− x) + q(x)

e valutandola in xp′(x) = q(x)

dove q(x) ed i suoi coefficienti sono quelli che si ottengono applicandol’algoritmo di Ruffini-Horner per valutare p(x) in x.

Esempio 1.2 Sia dato

p(x) = 1 + x− 2x2 + 3x4

e si voglia calcolare p(2) e p′(2)

3 0 −2 1 12 6 12 20 42

3 6 10 21 43 ≡ p(2)

2 6 24 683 12 34 89 ≡ p′(2)

Il calcolo di p(x) e p′(x) puo essere organizzato nel seguente modo:

p1:=0

p:=a[n]

per k=n-1,...,0

p1:=p+x_bar*p1

p:=a[k]+x_bar*p

p(x_bar):=p

p’(x_bar):=p1

Analogamente si possono calcolare anche le derivate di ordine superio-re. Derivando l’espressione ottenuta per p′(x) si ottiene:

p′′(x) = q′′(x)(x− x) + 2q′(x)

e valutando questa espressione in x

p′′(x) = 2q′(x).

Allora per ottenere p′′(x) e sufficiente calcolare q′(x) che possiamo otteneredallo schema di Horner utilizzato per valutare p′(x) in x.In generale se si indica la relazione fra polinomi dividendo, ed il suopolinomio quoziente, con un indice, cioe

pj(x) = pj+1(x− x) + pj(x) j = 0, . . . , n− 1

allora

1.2. VALUTAZIONE ED ERRORE INERENTE 7

• p(x) ≡ p0(x)

• p(j)(x) = j!pj(x)

• i coefficienti di pj sono i valori determinati applicando il metodo diHorner al polinomio pj−1.

1.2 Valutazione di un Polinomio: Errore Ine-

rente

Assegnato un polinomio ed un punto in cui valutarlo, l’errore inerentemisura a piccole variazioni sui coefficienti, come varia, in senso relati-vo, il valore del polinomio. Piu piccola e la variazione e migliore e ilcondizionamento del problema.

Esempio 1.3 Sia assegnato il polinomio

p(x) = a0 + a1x = 100− x

e lo si voglia valutare in punti x ∈ [100, 101]. Si perturba il coefficiente a1

dell’1% 1, per cui il polinomio perturbato risulta:

p(x) = 100− (1−1

100)x = 100−

99

100x

Valutandoli in x = 101 si ha:

p(101) = −1 p(101) = 0.01

commettendo un errore relativo dato da:∣

∣

∣

∣

∣

p(101)− p(101)

p(101)

∣

∣

∣

∣

∣

= 1011

100.

Risulta quindi che una perturbazione iniziale dell’1% porta una variazionesul risultato finale del 101%.La causa di questa amplificazione dell’errore e evidente nella Fig. 1.2.In pratica il coefficiente a1 rappresenta l’inclinazione della retta la quale,anche se alterata minimamente, comporta grossi errori per punti lontanidall’origine.Lo stesso comportamento si ottiene perturbando entrambi i coefficienti ovalutando in altri punti dell’intervallo.


Figura 1.2: Amplificazione dell’errore

In generale il problema di valutare un polinomio dato nella base deimonomi in punti appartenenti ad un intervallo [a, b] con a e b grandi ea/b ≃ 1 non risulta un problema ben condizionato, ossia si possono averegrossi errori inerenti.Procediamo, per questo esempio, all’analisi dell’errore inerente mediantela nota stima:

ǫin ≃n∑

i=1

ciǫi (1.4)

dove

ci =xi

f(x)

∂f(x)

∂xi

(1.5)

con ǫi = xi−xi

xi

, e |ǫi| < u, per i = 1, 2, . . . , n.Nel caso specifico di p(x) = a0 + a1x sara:

ǫin ≃ c1ǫ1 + c2ǫ2 + c3ǫ3

=a0

a0 + a1xǫ1 +

a1x

a0 + a1xǫ2 +

xa1

a0 + a1xǫ3

e nel caso a0 = 100, a1 = −1 e x ∈ [100, 101] con, per esempio, x = 101 siha:

=100

−1ǫ1 +

−101

−1ǫ2 +

−101

−1ǫ3

1cioe | a1−a1

a1

| = 1100 ; segue che a1 = a1 ±

1100a1 = (1± 1

100 )a1

1.2. VALUTAZIONE ED ERRORE INERENTE 9

da qui si vede come una piccola perturbazione su uno dei coefficienti, peresempio a1 dell’1% (ǫ1 = 0, ǫ2 = 1/100, ǫ3 = 0) porti ǫin ad assumere ilvalore −101 1

100e cioe 101 volte maggiore di quello iniziale.

Dall’analisi effettuata si evince che per qualunque punto x in [100, 101]questo comportamento sara inevitabile in quanto una lieve modifica deicoefficienti (ǫ1 o ǫ2 6= 0) viene grandemente amplificata (coefficienti c1 ec2). Si osservi inoltre che anche per ǫ1 = ǫ2 = 0 e ǫ3 6= 0 si avra un erroreinerente grande in questo intervallo di valutazione.

Generalizzando l’analisi fatta in questo esempio alla valutazione di ungenerico polinomio p(x) nella base monomiale, l’errore inerente si puopresentare in due componenti:

ǫin1 =n∑

i=0

ai

p(x)

∂p(x)

∂ai

ǫi =n∑

i=0

aixi

p(x)ǫi con ǫi =

ai − ai

ai

e

ǫin2 =x

p(x)

∂p(x)

∂xǫx =

xp′(x)

p(x)ǫx con ǫx =

x− x

x.

Si puo quindi desumere che:

• ǫin1 dipende da p(x) e dai valori aixi. Questi termini sono in parte

controllabili riformulando il problema (la sua rappresentazione);

• ǫin2 dipende unicamente dal polinomio in esame p(x) e dalla suaderivata, per cui e un errore intrinseco al valore stesso di p(x) e nondipende dalla sua rappresentazione.In particolare si ha un grande errore relativo per p(x)→ 0 o perp′(x)→∞.

Nel seguito vogliamo esaminare se e possibile agire sull’espressione di p(x)al fine di ridurre l’errore inerente. Si introduce il concetto di condizio-namento di una base polinomiale per indicare quale rappresentazione po-linomiale sia affetta da minor errore inerente allorquando se ne usi unasua combinazione lineare su un intervallo per risolvere un problema chepuo andare dalla valutazione alla determinazione dei suoi zeri. Si osserviche, come introdotto anche nell’esempio, lavorando numericamente restasempre definito l’intervallo di lavoro o di interesse e questa informazionepuo essere sfruttata per individuare una base meglio condizionata.


1.3 Polinomi nella base di Bernstein

Con polinomio in forma di Bernstein si intende una espressione del tipo

p(x) =n∑

i=0

biBi,n(x) con x ∈ [a, b]

dove {Bi,n(x)} sono le funzioni base di Bernstein definite da

Bi,n(x) =

(

n

i

)

(b− x)n−i(x− a)i

(b− a)ni = 0, . . . , n

e b0, b1, . . . , bn ∈ IR sono i coefficienti nella base di Berstein.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figura 1.3: Polinomi base di Bernstein di grado 3.

Esempio 1.4 Riprendiamo l’esempio del polinomio precedente, ma nellabase di Bernstein; sara:

B0,1(x) = 101− x B1,1(x) = x− 100

da cuip(x) = 100− x = b0B0,1(x) + b1B1,1

= 0(101− x)− 1(x− 100) ossia b0 = 0 b1 = −1.

Rieseguendo l’analisi sull’errore inerente si ha:

ǫin = c1ǫ1 + c2ǫ2 + c3ǫ3 =

b0

p(x)(101− x)ǫ1 +

b1

p(x)(x− 100)ǫ2 +

xp′(x)

p(x)ǫ3

1.3. POLINOMI DI BERNSTEIN 11

e perturbando solo b1 (cioe ǫ1 = ǫ3 = 0, ma ǫ2 6= 0), valutando in unqualunque punto in [100, 101] il coefficiente moltiplicatore di ǫ2 sara almassimo dell’ordine dell’unita.Procediamo numericamente; si perturba il coefficiente b1 dell’1%, per cuiil polinomio perturbato risulta:

p(x) = −(1−1

100)(x− 100) = −

99

100(x− 100)

Valutandoli in x = 101 si ha:

p(101) = −1 p(101) = −0.99


∣

∣

∣

∣

p(101)− p(101)

p(101)

∣

∣

∣

∣

∣

=1

100

Risulta quindi che una perturbazione iniziale dell’1% porta ad una varia-zione sul risultato finale della stessa entita.Come gia fatto notare anche nell’analisi precedente, se ǫ3 6= 0, xp′(x)

p(x)po-

trebbe essere grande indipendentemente dalla base di rappresentazione.Vediamo un esempio numerico.Perturbiamo x dell’1%; se x = 101 sara x = 99

100101 = 99.99; allora:

p(101) = −1 p(99.99) = 0.01


∣

∣

∣

∣

p(101)− p(99.99)

p(101)

∣

∣

∣

∣

∣

= 1011

100;

si puo fare qualcosa?

Si ricorda che i polinomi godono della proprieta di essere invarianti pertraslazione e scala dell’intervallo di definizione o cambio di variabile. Ciosignifica che dato un polinomio p(x) in un intervallo [a, b] e un intervallotraslato e scalato di questo, per esempio [0, 1], e possibile definire unaapplicazione (mapping) dal primo al secondo e viceversa

x ∈ [a, b]→ t ∈ [0, 1]

t =x− a

b− ao viceversa x = a + t(b− a) (1.6)

ed effettuando un cambio di variabile determinare un polinomio q(t) cheper ogni t ∈ [0, 1] assume gli stessi valori di p(x) per la x corrispondente


secondo il mapping definito. Un problema nell’applicare tale cambio divariabile consiste nel dover determinare numericamente i coefficienti delpolinomio q(t) rischiando di commettere degli errori di approssimazione.I polinomi di Bernstein godono della proprieta che nel cambio di variabilei coefficienti non cambiano. Vediamolo nel dettaglio.

Bi,n(x) = Bi,n(a+t(b−a)) =

(

n

i

)

[b− a− t(b− a)]n−i[a + t(b− a)− a]i

(b− a)n=

=

(

n

i

)

[(b− a)(1− t)]n−i[t(b− a)]i

(b− a)n=

(

n

i

)

(1− t)n−iti = Bi,n(t) (1.7)

Riassumendo, avremo

p(x) =n∑

i=0

biBi,n (x) x ∈ [a, b]

ed effettuando il cambio di variabile otteniamo

p(t) =n∑

i=0

biBi,n (t) t ∈ [0, 1]

che dice che: un cambio di variabile, per un polinomio nella base di Bern-stein, non comporta alcun errore o costo computazionale sui coefficientiperche restano i medesimi.

Tornando all’errore inerente nella valutazione ed in particolare al no-stro esempio numerico si ha:

p(x) = b0(101− x) + b1(x− 100) x ∈ [100, 101]

p(t) = b0(1− t) + b1(t− 0) t ∈ [0, 1]

e la stima dell’errore inerente

ǫin =b0

p(t)(1− t)ǫ1 +

b1

p(t)(t− 0)ǫ2 +

tp′(t)

p(t)ǫ3

dove i primi due termini restano uguali in valore, mentre il terzo purrestando p(t) lo stesso, potenzialmente puo avere tp′(t) inferiore; questo esicuramente vero per t che assume valori in [0, 1] e se si richiama che valela relazione p′(t) = (b − a)p′(x), verra ridotto anche p′(t) se (b − a) < 1.Tornando all’esempio numerico si ha:

p(t) = −t t ∈ [0, 1]


e poiche in questo caso b− a = 1 vale p′(t) = p′(x). Peturbiamo t dell’1%,allora t = 99

100t; se t = 1 sara t = 0.99, allora

p(1) = −1 p(0.99) = −0.99


∣

∣

∣

∣

p(1)− p(0.99)

p(1)

∣

∣

∣

∣

∣

=1

100.

Se b − a > 1, si puo suddividere l’intervallo di interesse in piu intervallitutti di ampiezza minore di 1 (vedi sezione 1.3.5).Si noti che in [0, 1] i polinomi di Bernstein sono piu semplici da calcolare(vedi la relazione 1.7 che li definisce) e la valutazione del polinomio risultain genere piu stabile dal punto di vista numerico.Riassumendo, se si vuole valutare p(x), si determina t := ¯x−a

(b−a), si calcola

p(t) con p(t) polinomio di Bernstein in [0, 1], e sara p(x) ≡ p(t).Da ora in poi useremo sempre i polinomi di Bernstein in [0, 1].

1.3.1 Errore algoritmico: un esempio famoso

Della valutazione numerica di un polinomio abbiamo attentamente ana-lizzato le cause che possono portare ad un errore inerente grande. Anchese solo con un esempio si vuole evidenziare l’errore algoritmico dovutoa propagazione di errori e mostrare che l’arrotondamento puo rovinareuna computazione ordinaria. A tal fine considereremo un polinomio nel-la base monomiale con coefficienti esattamente rappresentabili (errore dirappresentazione nullo) e utilizzeremo il metodo di Horner per valutaretale polinomio in punti di un intervallo che siano numeri finiti (errore dirappresentazione nullo sui punti di valutazione). Grazie poi al fatto chetale polinomio ammette una rappresentazione esatta anche nella base diBernstein, potremo utilizzare un metodo alternativo per la sua valutazionee confrontare le risposte dei due metodi.

Esempio 1.5 Si consideri il polinomio

p(x) = 1− 6x + 15x2 − 20x3 + 15x4 − 6x5 + x6

e lo si valuti, usando la precisione basic double, nell’intervallo [a, b] =[0.99609375, 1.00390625] di estremi numeri finiti rappresentabili in base 2con meno di 8 cifre, usando un passo h = 0.000244140625 per un totale


di 33 punti rappresentabili in base 2 tutti al piu con 12 cifre. Lo stessopolinomio ammette la seguente rappresentazione nella base di Bernstein:

p(x) =

{

(1− x)6 x ∈ [0, 1](x− 1)6 x ∈ [1, 2]

La Fig. 1.4 mostra i grafici delle valutazioni effettuate con la rappresen-tazione monomiale e quella di Bernstein ed evidenzia errori algoritmicinell’utilizzo del metodo di Horner del polinomio in forma monomiale. La

0.995 0.996 0.997 0.998 0.999 1 1.001 1.002 1.003 1.004 1.005−0.5

0

0.5

1

1.5

2

2.5

3

3.5

4x 10

−15

Figura 1.4: Grafici delle valutazioni effettuate in 33 punti.

Fig.1.5 mostra i grafici di analoghe valutazioni, ma in 201 punti equidi-stanti dell’intervallo [a, b] = [0.995, 1.005]; ancora maggiormente si notal’instabilita della valutazione in oggetto.

Osservazione 1.1 Le basi di rappresentazione polinomiale sono infinitee numerose sono quelle di interesse numerico, nel senso che a seconda delproblema, si puo individuare la base piu idonea per trattarlo; fra questericordiamo la base delle potenze traslate con un centro, la forma di Newton(o base delle potenze traslate con n centri), la base dei polinomi cardinalio di Lagrange ed altre ancora. Forse la piu eclettica per le sue svariatebuone proprieta e ottime qualita numeriche e proprio la base di Bernstein,che abbiamo deciso di sponsorizzare in questa trattazione.


0.995 0.996 0.997 0.998 0.999 1 1.001 1.002 1.003 1.004 1.005−4

−2

0

2

4

6

8

10

12

14

16x 10

−15

Figura 1.5: Grafici delle valutazioni effettuate in 201 punti.

1.3.2 Proprieta dei polinomi di Bernstein

I polinomi base di Bernstein Bi,n(t) godono di alcune interessanti pro-prieta:

• Bi,n(t) ≥ 0 i = 0, . . . , n ∀t ∈ [0, 1];

• Partizione dell’unita:n∑

i=0

Bi,n(t) = 1 ∀t ∈ [0, 1]

• Per le proprieta precedenti, p(t) espresso nella base di Bernstein euna combinazione convessa dei bi, da cui segue

mini{bi} ≤ p(t) ≤ max

i{bi} ∀t ∈ [0, 1].

Si richiama che una combinazione lineare si dice affine se la sommadei coefficienti e 1 (

∑

Bi,n(t) = 1), se poi questi sono non negativi(0 ≤ Bi,n(t) ≤ 1) allora si dice combinazione convessa.

• Variazione di segno dei coefficienti:un polinomio espresso nella base di Bernstein ha un numero di varia-zioni di segno minore o uguale al numero di variazione di segno delvettore dei suoi coefficienti (valori o coefficienti nulli non vengonoconsiderati).


1.3.3 Formula ricorrente

I polinomi base di Bernstein di grado n, sono legati ai polinomi base diBernstein di grado n− 1, dalla seguente formula ricorrente

Bi,n(t) = tBi−1,n−1(t) + (1− t)Bi,n−1(t)

con B0,0(t) = 1 e Bi,n(t) = 0 ∀i /∈ [0, n].

Viene riportato il codice Matlab che implementa la formula ricorrente,sui polinomi base di Bernstein in [0,1], per la valutazione in un punto o inun array di punti.

function bs=bernst(n,t)

% calcola i polinomi di Bernstein in [0,1] in un punto t;

% se t e’ un vettore di punti torna una matrice di valori;

% n --> grado del polinomio

% x --> valore di un punto o vettore di punti

% bs <-- matrice dei polinomi di bernstein nei punti

m=length(t);

np1=n+1;

bs=zeros(m,np1);

for ii=1:m

l=np1;

bs(ii,l)=1.0;

d1=t(ii);

d2=1.0-t(ii);

for i=1:n

temp=0.0;

for j=l:np1

beta=bs(ii,j);

bs(ii,j-1)=d2.*beta+temp;

temp=d1.*beta;

end

bs(ii,np1)=temp;

l=l-1;

end

end

Osservazione 1.2 Mediante tale algoritmo la valutazione dei polinomibase di Bernstein in un punto costa n(n+1)

2addizioni/sottrazioni ed n(n+1)

moltiplicazioni.


Esempio 1.6 Polinomi base di Bernstein di grado 3 ottenuti medianteformula ricorrente:

B0,0(t) = 1B0,1(t) = (1− t) B1,1(t) = tB0,2(t) = (1− t)2 B1,2(t) = 2t(1− t) B2,2(t) = t2

B0,3(t) = (1− t)3 B1,3(t) = 3t(1− t)2 B2,3(t) = 3t2(1− t) B3,3(t) = t3

1.3.4 Valutazione via algoritmo di de Casteljau

Per valutare un polinomio nella base di Bernstein si puo utilizzare laformula ricorrente sulle funzioni base per valutarle nel punto desiderato,quindi effettuare la combinazione convessa

p(xv) =n∑

i=0

biBi,n(xv).

Alternativamente si puo usare l’algoritmo di de Casteljau sui coefficientidel polinomio, che deriva dall’applicare ripetutamente la formula ricorren-te vista.

p(t) =n∑

i=0

biBi,n(t)def=

n∑

i=0

bitBi−1,n−1(t) +n∑

i=0

bi(1− t)Bi,n−1(t) =

poiche B−1,n−1 = Bn,n−1 = 0

=n−1∑

i=0

bi+1tBi,n−1(t) +n−1∑

i=0

bi(1− t)Bi,n−1(t) =

raccogliendo

=n−1∑

i=0

[bi+1t + bi(1− t)]Bi,n−1(t) =

=n−1∑

i=0

b[1]i Bi,n−1(t) =

riapplicando la formula ricorrente piu volte, alla fine si ottiene:

0∑

i=0

b[n]0 B0,0(t) = b

[n]0 .

Si ha quindi il seguente algoritmo:

b[k]i = tb

[k−1]i+1 + (1− t)b

[k−1]i (1.8)


con k = 1, . . . , n, i = 0, . . . , n − k e p(t) := b[n]0 ; in modo esplicito si ha il

seguente schema triangolare:

b[0]0 b

[0]1 b

[0]2 . . . b

[0]n−1 b[0]

n

b[1]0 b

[1]1 b

[1]2 . . . b

[1]n−1

. . . . . .

b[n−2]0 b

[n−2]1 b

[n−2]2

b[n−1]0 b

[n−1]1

b[n]0

(1.9)

Osservazione 1.3 Mediante tale algoritmo la valutazione di un polino-mio nella base di Bernstein costa n(n+1)

2addizioni/sottrazioni ed n(n + 1)

moltiplicazioni.

function y=decast(c,t)

% calcola il valore di un polinomio nella base di Bernstein

% definito in [0,1] dai coefficienti c nei punti t mediante

% l’algoritmo di de Casteljau;

% c --> coefficienti del polinomio

% t --> vettore dei punti

% y <-- vettore dei valori del polinomio nei punti

np1=length(c);

n=np1-1;

m=length(t);

for k=1:m

w=c;

d1=t(k);

d2=1.0-t(k);

for j=1:n

for i=i:np1-j

w(i)=d1.*w(i+1)+d2.*w(i);

end

end

y(k)=w(1);

end

Osservazione 1.4 Il valore di un polinomio di Bernstein agli estremi ebanalmente dato dai suoi coefficienti b0 e bn, cioe:

p(0) := b0 e p(1) := bn.


1.3.5 Suddivisione

Una utile applicazione dell’algoritmo di de Casteljau e quella di poter es-sere utilizzato per suddividere un polinomio p(t) in un punto tc, ottenendoi suoi coefficienti nelle basi di Bernstein per gli intervalli [0, tc] e [tc, 1]; talicoefficienti risultano essere rispettivamente:

b[0]0 , b

[1]0 , . . . , b

[n]0 e b

[n]0 , b

[n−1]1 , . . . , b[0]

n (1.10)

che corrispondono ai lati verticale e diagonale dello schema triangolare1.9. Si noti che, grazie all’invarianza per traslazione e scala vista, questirestano gli stessi anche quando i due intervalli vengono mappati in [0, 1]mediante le trasformazioni

t

tc≡ u ∈ [0, 1] e

t− tc1− tc

≡ v ∈ [0, 1].

Per dimostrare quanto asserito, notiamo che le quantita b[k]i generate me-

diante la formula ricorrente 1.8 possono essere espresse in termini deicoefficienti iniziali bi nella base di Bernstein nella forma:

b[k]i =

k∑

j=0

bi+jBj,k(tc) i = 0, . . . , k

che e facilmente verificata per come sono stati costruiti. I coefficienti 1.10possono allora essere scritti come:

b[k]0 =

k∑

j=0

bjBj,k(tc) e b[n−k]k =

n−k∑

j=0

bk+jBj,n−k(tc) k = 0, . . . , n. (1.11)

Si noti che la suddivisione descritta altro non e che un cambio di base daBernstein in [0, 1] a Bernstein in [0, tc] e [tc, 1] rispettivamente. A questoproposito si richiama il seguente risultato.

Teorema 1.4 Sia {Bj,n(t)} la base di Bernstein in [0, 1], allora le basi diBernstein in [0, tc] e [tc, 1], che indicheremo rispettivamente con {Bj,n(t)}e {Bj,n(t)}, sono date da:

[B0,n, B1,n, . . . , Bn,n] = [B0,n, B1,n, . . . , Bn,n]Γ

[B0,n, B1,n, . . . , Bn,n] = [B0,n, B1,n, . . . , Bn,n]Λ

con

Γ =

B0,0(tc) 0 . . . 0B0,1(tc) B1,1(tc) . . . 0. . .B0,n(tc) B1,n(tc) . . . Bn,n(tc)


e

Λ =

B0,n(tc) B1,n(tc) . . . Bn,n(tc). . .0 . . . B0,1(tc) B1,1(tc)0 . . . 0 B0,0(tc)

dim.Se sostituiamo t = tcu, 1 − t = (1 − u) + (1 − tc)u e t = tc(1 − v) + v,1 − t = (1 − tc)(1 − v) nella definizione 1.7 della funzione base Bj,n(t)t ∈ [0, 1], si ottiene il seguente sviluppo in termini delle basi Bi,n(u) eBi,n(v) su u ∈ [0, 1] e v ∈ [0, 1] rispettivamente.

Bj,n(t) =

(

n

j

)

(1− t)n−jtj =

=

(

n

j

)

[(1− u) + (1− tc)u]n−j[tcu]j

e per lo sviluppo della potenza di un binomio

=

(

n

j

)

n−j∑

k=0

(

n− j

k

)

(1− u)n−j−k(1− tc)kuktjcu

j

shiftando gli indici della sommatoria si ha

=

(

n

j

)

n∑

k=j

(

n− j

k − j

)

(1− u)n−k(1− tc)k−juk−jtjcu

j

e osservando che(

n

j

)(

n− j

k − j

)

≡

(

n

k

)(

k

j

)

si ha

=n∑

k=j

(

k

j

)

(1− tc)k−jtjc

(

n

k

)

(1− u)n−kuk

=n∑

k=j

Bj,k(tc)Bk,n(u)

che dimostra il primo cambio di base; per il secondo si procede analoga-mente.

E ora banale osservare che quanto si ottiene con l’algoritmo di de Ca-steljau applicato in un punto tc e che si e riscritto nella forma 1.11, corri-sponde proprio ai coefficienti della rappresentazione dello stesso polinomio,


ma nelle basi di Bernstein in [0, tc] e [tc, 1]; infatti da

p(t) =n∑

i=0

biBi,n(t) = [B0,n(t), B1,n(t), . . . , Bn,n(t)]

b0

b1

· · ·bn

applicando il cambio di base, per esempio in [0, tc] si ha

= [B0,n(u), B1,n(u), . . . , Bn,n(u)]Γ

b0

b1

· · ·bn

dove

Γ

b0

b1

· · ·bn

equivale proprio alla rappresentazione della prima delle 1.11; analogamen-te per la seconda.

Si noti che ogni suddivisione di un polinomio raddoppia il numero dicoefficienti di Bernstein richiesti per rappresentarlo sull’intero intervallooriginale. Applicando un processo iterativo di suddivisione uniforme ivalori dei coefficienti su ciascun sottointervallo convergono ai valori delpolinomio. Questa osservazione e alla base della procedura per generareuna approssimazione lineare a tratti di un polinomio di Bernstein che asua volta e alla base di numerose applicazioni.

Si noti inoltre che una suddivisione, dal punto di vista del costo com-putazionale, altro non e che una valutazione in cui si memorizzano al-cuni coefficienti intermedi. L’idea di valutare/suddividere in un puntodell’intervallo, quindi utilizzare i coefficienti della suddivisione per valu-tare/suddividere ulteriormente, fa sı che il condizionamento per ogni va-lutazione succesiva migliori. Questo e dovuto semplicemente al fatto cheogni sottointervallo, una volta operata una suddivisione, viene rimappatoin [0, 1] per la successiva valutazione/suddivisione.

1.3.6 Derivata

La derivata di un polinomio p(t) nella base di Bernstein sara espressa da:

p′(t) =n∑

i=0

biB′

i,n(t)


dove

B′

i,n(t) = n(Bi−1,n−1(t)−Bi,n−1(t)).

Sostituendo si ottiene:

p′(t) = nn∑

i=0

bi[Bi−1,n−1(t)−Bi,n−1(t)] =

= nn∑

i=0

biBi−1,n−1(t)− nn∑

i=0

biBi,n−1(t) =

= nn−1∑

i=0

bi+1Bi,n−1(t)− nn−1∑

i=0

biBi,n−1(t) =

=n−1∑

i=0

diBi,n−1(t)

con di = n(bi+1− bi) i = 0, . . . , n−1 che sono i coefficienti del polinomioderivata prima, di grado n− 1, nella base di Bernstein.

Osservazione 1.5 I valori della derivata di un polinomio di Bernsteinagli estremi sono banalmente dati dai coefficienti d0 e dn−1, cioe:

p′(0) := d0 e p′(1) := dn−1.

Osservazione 1.6 La derivata di un polinomio di Bernstein di grado ne definito su [a, b] e data da:

p′(x) =1

b− ap′(t)

con p′(t) la derivata di p(t) in [0, 1].

Osservazione 1.7 Nel caso la situazione richiedesse contemporaneamen-te il valore ed il valore della derivata prima in un punto t, si puo procederenel seguente modo: si applichi l’algoritmo di de Casteljau per valutare ilpolinomio ottenendo p(t) = b

[n]0 , quindi la derivata prima sara data da

p′(t) = n(b[n]0 − b

[n−1]0 ) ≡ n(b

[n−1]1 − b

[n]0 ); questo deriva dal fatto che l’al-

goritmo di de Casteljau nel valutare in t il polinomio ne fornisce unarappresentazione nella base di Bernstein negli intervalli [0, t] e [t, 1] e ne-gli estremi dell’intevallo di definizione il valore ed il valore della derivatasono noti in modo banale.


function [y,yd]=decast_der(c,t)

% calcola il valore ed il valore della derivata di un polinomio

% nella base di Bernstein definito in [0,1) dai coefficienti c

% nei punti t mediante l’algoritmo di de Casteljau;

% c --> coefficienti del polinomio

% t --> vettore dei punti

% y <-- vettore dei valori del polinomio nei punti

% yd <-- vettore dei valori di derivata del polinomio nei punti

np1=length(c);

n=np1-1;

m=length(t);

for k=1:m

w=c;

d1=t(k);

d2=1.0-t(k);

for j=1:n

for i=i:np1-j

w(i)=d1.*w(i+1)+d2.*w(i);

end

end

y(k)=w(1);

yd(k)=n(w(2)-w(1))/d2;

end

1.3.7 Antiderivata e integrazione

Sia p(t) =∑n−1

i=0 diBi,n−1(t) un polinomio di grado n−1 definito in [0, 1] nel-la base di Bernstein. Una sua antiderivata o primitiva sara un polinomiodi grado n che nella base di Bernstein puo essere espressa come

p−1(t) =n∑

i=0

biBi,n(t)

con

bi+1 :=di

n+ bi i = 0, . . . , n− 1

avendo fissato un valore per b0 (per esempio b0 := 0).La scelta arbitraria del valore iniziale b0 e data dal fatto che esistono piuprimitive tutte differenti a meno di una costante additiva.

Volendo procedere a determinare l’integrale definito di un polinomiop(t) si avra:

∫ 1

0p(t)dt = [p−1(t)]10 = bn − b0 = bn


se, come suggerito prima, si sceglie b0 = 0. In particolare sara poi:

bn =dn−1

n+ bn−1 ==

dn−1

n+

dn−2

n+ bn−2 = . . . =

1

n

n−1∑

i=0

di.

Osservazione 1.8 L’integrale definito di ogni polinomio di Bernstein digrado n in [0, 1] e la costante 1

n+1.

Infatti Bi,n(t) e un polinomio nella base di Bernstein di coefficienti

dj =

{

0 j 6= i1 j = i

cosı che∫ 1

0Bi,n(t)dt =

1

n + 1

n∑

j=0

dj =1

n + 1.

Osservazione 1.9 L’integrale definito di ogni polinomio base di Bern-stein di grado n in [a, b] e la costante b−a

n+1.

Questo e banalmente vero pensando di calcolare

∫ b

aBi,n(x)dx

effettuando un cambio di variabile da [a, b] a [0, 1]; infatti sara:

∫ b

aBi,n(x)dx = (b− a)

∫ 1

0Bi,n(t)dt =

b− a

n + 1

1.4 Un’applicazione: le curve di Bezier

Bezier, negli anni ′60, fu l’ideatore di uno dei primi sistemi CAD, UNI-SURF [Bez70], usato dalla casa automobilistica francese Renault. La teo-ria delle curve di Bezier, su cui era basato UNISURF, nasce dal rappre-sentare un polinomio nella base di Bernstein.Una curva piana di Bezier C(t) di grado n puo essere definita a partireda un insieme di punti di controllo Pi ∈ E2 (piano Euclideo) i =0, . . . , n ed e data da

C(t) =n∑

i=0

PiBi,n(t) t ∈ [0, 1].

Si consideri per esempio una curva di grado 3. La Fig. 1.6 a sinistramostra una tale curva e i punti di controllo associati P0, P1, P2 e P3. Si

1.4. UN’APPLICAZIONE: LE CURVE DI BEZIER 25

Figura 1.6: Esempio di curva di Bezier; caso n = 3 a sinistra; modifica diforma mediante spostamento di un punto di controllo a destra.

puo vedere semplicemente che i punti estremi della curva coincidono conP0 e P3 e che i punti di controllo non sono punti della curva.Spostando i punti di controllo si influenza la forma della curva; la Fig. 1.6a destra mostra l’effetto di spostare il punto P1 in P ′

1.E’ facile pensare ad un software interattivo che permetta ad un utente dimuovere i punti di controllo e quindi deformare la curva in una manierache risulti subito intuitiva. Questa semplice possibilita offerta dai puntidi controllo e il fondamento dell’importanza e utilita pratica delle curvedi Bezier.

Le curve di Bezier sono oggi utilizzate in molti sistemi CAD 2D (Com-puter Aided Design), CAD 3D, sistemi di disegno, software e linguaggigrafici come Adobe Illustrator, Postscript, OpenGL, METAFONT, La-Tex, Xfig, gimp, inkbase e tanti altri. Le curve di Bezier sono poi unostandard per la rappresentazione di fonti di caratteri, di modelli CAD ingenere e sono alla base dello standard SVG (Scalable Vector Graphics)per la grafica vettoriale in alternativa alla grafica raster.

Capitolo 2

Interpolazione

Siano assegnate n+1 osservazioni in corispondenza di n+1 punti distinti,cioe siano assegnati (xi, yi) i = 0, . . . , n. Nel caso generale viene indivi-duato uno spazio a dimensione finita Φ di funzioni reali e una loro base dirappresentazione φ0, φ1, . . . , φn cosı che ogni elemento φ ∈ Φ sia definitoda una loro combinazione lineare, cioe

φ =n∑

i=0

aiφi(x).

Allora il problema di interpolare i dati assegnati con una funzione φ ∈ Φconsiste nel determinare la funzione φ∗ (od anche i coefficienti a∗

0, a∗1, . . . , a

∗n

che la definiscono) tale che

φ∗(xi) = yi i = 0, . . . n.

Spesso, in pratica, questo problema viene applicato per determinareuna approssimazione su un intervallo di una data funzione; in questo casoviene assegnata una funzione continua f(x), x ∈ [a, b] e la si vuole ap-prossimare con una funzione, in genere piu semplice, φ ∈ Φ, ma in modoabbastanza fedele cosı che la funzione errore

R(x) = |f(x)− φ(x)| < ǫ x ∈ [a, b]

con ǫ una costante sufficientemente piccola legata all’applicazione.Nel caso generale qui descritto, la scelta dello spazio Φ e solitamente

dettata dal fatto che esista sempre una soluzione unica al problema, dallaregolarita desiderata per l’interpolante e dalla particolarita dei dati.

In queste dispense si trattera solamente il caso dell’interpolazione po-linomiale che garantisce sempre l’esistenza e l’unicita della soluzione, unabuona regolarita ed una buona ricostruzione in casi particolari.

27

28 CAPITOLO 2. INTERPOLAZIONE

2.1 Forma Monomiale

Teorema 2.1 Siano dati n + 1 punti (xi, yi), i = 0, . . . , n con xi distinti,allora esiste ed e unico il polinomio p ∈ Pn che verifica le condizioni

p(xi) = yi i = 0, . . . , n

dim. Si consideri p(x) nella forma monomiale

p(x) = a0 + a1x + . . . + anxn

e si impongano le n + 1 condizioni di interpolazione

a0 + a1x0 + a2x20 . . . + anx

n0 = y0

a0 + a1x1 + a2x21 . . . + anx

n1 = y1

. . ....

a0 + a1xn + a2x2n . . . + anx

nn = yn

Si tratta di un sistema lineare di n+1 equazioni in n+1 incognite; se talesistema ammette una ed una sola soluzione, allora tale sara il polinomiodi grado ≤ n. In forma matriciale il sistema si presentera come:

V a = y (2.1)

con

V =

1 x0 x20 . . . xn

0

1 x1 x21 . . . xn

1

. . .1 xn x2

n . . . xnn

a =

a0

a1...an

y =

y0

y1...yn

dove V e nota come matrice di Vandermonde; il determinante di talematrice ha la seguente espressione analitica:

detV =n∏

i,j=0,j>i

(xj − xi)

che nelle ipotesi che i punti xi siano distinti risulta sempre non nullo.Segue che il sistema lineare ha un’unica soluzione e quindi p(x) esiste ede unico.

La dimostrazione del Teorema fornisce anche un metodo per proce-dere alla determinazione del polinomio interpolante risolvendo il sistemalineare 2.1; purtroppo questo non risulta numericamente stabile a cau-sa del mal condizionamento del problema di determinare la soluzione delsistema lineare con matrice dei coefficienti la matrice di Vandermonde.

2.2. FORMA DI BERNSTEIN 29

2.2 Forma di Bernstein

Nella sezione precedente si e dimostrata l’esistenza ed unicita del polino-mio interpolante nella base monomiale, ma il polinomio interpolante puoessere rappresentato in una qualunque altra base. Procediamo quindi allasua determinazione mediante soluzione di un sistema lineare, ma a partiredalla base di Bernstein. Siano dati n + 1 punti (xi, yi), i = 0, . . . , n conxi distinti in [a, b] e si vuole determinare p(x) nella base di Bernstein taleche

p(xi) = yi i = 0, . . . , n

Siano (ti, yi), i = 0, . . . , n i punti traslati e scalati in [0, 1] mediante la 1.6e sia

p(t) =n∑

i=0

biBi,n(t) con t ∈ [0, 1]

la forma polinomiale che vogliamo interpoli i valori assegnati. Imponendole condizioni di interpolazione si ottiene il seguente sistema lineare:

Bb = y

con

B =

B0,n(t0) B1,n(t0) . . . Bn,n(t0)B0,n(t1) B1,n(t1) . . . Bn,n(t1). . .B0,n(tn) B1,n(tn) . . . Bn,n(tn)

b =

b0

b1...bn

y =

y0

y1...yn

Il vettore b, soluzione del sistema lineare sopra dato fornira sia il polinomiointerpolante in [0, 1] che in [a, b]. La soluzione di un sistema lineare n× ncosta O(n3/3) operazioni floating point.

Si noti che la matrice dei polinomi base di Bernstein nei punti risultameglio condizionata che non la matrice di Vandermonde, come si puo osser-vare nella Tab. 2.2 in cui vengono mostrati gli indici di condizionamentodelle matrici calcolate su punti equidistanti nell’intervallo di definizionemediante polinomi di Bernstein e base delle potenze.

2.3 Forma di Newton

Come accennato nell’introduzione a questo capitolo, a seconda del proble-ma che si vuole risolvere, si puo determinare una base piu idonea delle altre


grado n Bernst. Potenze3 2.27 1.41E024 4.52 9.04E025 9.65 3.78E036 21.60 2.50E047 49.98 1.30E058 118.38 8.10E059 285.34 5.08E0610 697.14 3.04E07

Tabella 2.1: Indici di condizionamento delle matrici ottenute valutando ipolinomi di Bernstein e la base delle potenze in punti equidistanti.

sia dal punto di vista del condizionamento, ma anche che permetta di usa-re un metodo piu efficiente oltre che stabile. Nel caso dell’interpolazionepolinomiale sicuramente la base piu idonea e quella di Newton.

Siano dati i punti xi, i = 0, . . . , n distinti, allora le funzioni

1, (x− x0), (x− x0)(x− x1), . . . , (x− x0)(x− x1) · · · (x− xn−1)

formano una base polinomiale per lo spazio Pn detta base di Newton.Assegnati (xi, yi), i = 0, . . . , n, si voglia determinare il polinomio interpo-lante nella base di Newton, cioe si voglia determinare il polinomio p(x)nella forma

p(x) = c0+c1(x−x0)+c2(x−x0)(x−x1)+. . . cn(x−x0)(x−x1) · · · (x−xn−1).

Imponendo le condizioni di interpolazione si ottiene il sistema lineare:

Nc = y

con

N =

1 0 0 . . . 01 (x1 − x0) 0 . . . 01 (x2 − x0) (x2 − x0)(x2 − x1) . . . 0...1 (xn − x0) (xn − x0)(xn − x1) . . . (xn − x0) · · · (xn − xn−1)

c =

c0

c1...cn

y =

y0

y1...yn

2.4. FORMA DI LAGRANGE 31

che risulta facilmente risolubile per sostituzione in avanti con una com-plessita di O(n2). Tale soluzione e ottimale da ogni punto di vista, inoltree applicabile una variante dell’algoritmo di Horner per la sua valutazione.Purtroppo se il polinomio cosı determinato deve essere sottoposto ad ul-teriori elaborazioni come per esempio a derivazione o integrazione, alloratale forma non e piu raccomandabile in quanto non sono noti algoritmiefficienti per procedere.

2.4 Forma di Lagrange

Con forma di Lagrange si intende un polinomio nella base delle funzionicardinali o polinomi elementari di Lagrange. Siano dati i punti xi, i =0, . . . , n distinti, allora le funzioni elementari di Lagrange sono i polinomiLi,n(x) di grado n che soddisfano le seguenti condizioni:

Li,n(xj) =

{

1 se i = j0 se i 6= j

Questa base risulta la piu semplice in cui risolvere un problema di inter-polazione; assegnati (xi, yi), i = 0, . . . , n, il problema e risolto banalmenteda

p(x) =n∑

i=0

yiLi,n(x)

infatti per le condizioni che definiscono i polinomi Li,n(x), banalmente inogni punto xi l’espressione polinomiale data varra yi.Dal punto di vista computazionale determinare tale polinomio non costanulla essendo i coefficienti del polinomio in tale base proprio i valori yi as-segnati. Gia la valutazione del polinomio interpolante comporta la deter-minazione, ma soprattutto la valutazione dei polinomi Li,n(x). Per il Teo-rema 1.2 i polinomi Li,n(x) avendo come radici i punti xj con j = 0, . . . , ne j 6= i saranno della forma

Li,n(x) = ci(x− x0)(x− x1) · · · (x− xi−1)(x− xi+1) · · · (x− xn)

con ci una costante da determinare in modo che Li,n(xi) = 1; allora

ci =1

(xi − x0)(xi − x1) · · · (xi − xi−1)(xi − xi+1) · · · (xi − xn)

da cui in definitiva

Li,n(x) =

∏nj=0j 6=i(x− xj)

∏nj=0j 6=i(xi − xj)

i = 0, . . . , n


La valutazione di questo polinomio interpolante puo essere resa competiti-va con il polinomio di interpolazione nella forma di Newton, ma come quel-lo risulta difficoltosa la valutazione della derivata. Come si vedra nel capi-tolo dedicato all’integrazione numerica di funzioni, tale rappresentazionerisultera molto importante.

2.5 Errore nell’interpolazione polinomiale

Assegnate le n+1 osservazioni yi, i = 0, . . . , n in corrispondenza dei puntidistinti xi, si e visto come costruire il polinomio interpolante di gradominimo p(x). Se i valori yi altro non sono che i valori di una funzione f(x)nei punti xi, cioe yi ≡ f(xi), i = 0, . . . , n con f definita in [a, b], ha sensochiedersi quanto sia grande l’errore di interpolazione

R(x) = f(x)− p(x)

che si commette in un punto x ∈ [a, b] diverso dai punti di interpolazionexi. Per dare una risposta a tale domanda e necessario fare alcune ipotesi diregolarita sulla funzione f(x); infatti se f(x) non e almeno continua, l’er-rore fra i punti di interpolazione puo essere qualunque, anche grandissimo.Si dimostra il seguente teorema.

Teorema 2.2 Sia f(x) ∈ C(n+1)[a,b] , e sia x un punto qualsiasi in [a, b].

Allora esiste un punto ξ (dipendente da x) interno ad [a, b] per cui

f(x) = p(x) +w(x)

(n + 1)!f (n+1)(ξ)

con w(x) = (x− x0)(x− x1) · · · (x− xn).

Corollario 2.1 Posto M(n+1) = maxx∈[a,b] |f(n+1)(x)|, un limite superiore

all’errore di interpolazione R(x) = f(x)− p(x) e dato da

|R(x)| ≤|w(x)|

(n + 1)!M(n+1)

Esempio 2.1 Sia f(x) = ex; per ogni x ∈ [a, b] si ha M(n+1) = eb e perogni scelta dei punti xi, |w(x)| ≤ (b− a)(n+1), da cui

maxx∈[a,b]

|R(x)| ≤(b− a)(n+1)

(n + 1)!eb.

2.5. ERRORE NELL’INTERPOLAZIONE POLINOMIALE 33

In questo caso si ricava

limn→∞{max

x∈[a,b]|R(x)|} = 0,

cioe il polinomio interpolante p(x) converge a f(x) uniformemente su [a, b]quando il numero n dei punti di interpolazione tende all’infinito.

Osservazione 2.1 E opportuno rilevare che, anche se la funzione f(x) ∈C∞

[a,b], la successione pn(x) dei valori assunti dai polinomi di interpolazionedi grado n in un punto x ∈ [a, b] puo non convergere ad f(x), per n chetende all’infinito.

Esempio 2.2 Per il polinomio di interpolazione in n punti equidistantidella funzione di Runge

f(x) =1

1 + x2x ∈ [−5, 5]

si puo dimostrare che al crescere di n la successione dei polinomi di inter-polazione sui punti

xi =10

ni− 5

non converge puntualmente ad f(x) e che i corrispondenti resti diventanoin modulo arbitrariamente grandi in punti dell’intervallo [−5, 5].La Tab. 2.2 mostra l’errore massimo in valore assoluto fra la funzionedi Runge e i polinomi interpolanti. In particolare l’errore e grande vicinoagli estremi ed aumenta con n (vedi Fig.2.1).

Risultato 2.1 Se oltre a fare delle ipotesi sulla regolarita della funzione(almeno C1

[a,b]), si utilizzano opportune distribuzioni dei punti di interpo-lazione (come per esempio gli zeri del polinomio di Chebishev di gradon + 11), si puo dimostrare la convergenza del polnomio interpolante pn(x)alla f(x) per x ∈ [a, b] quando n → ∞ (si veda Fig.2.2 e la si confronticon Fig.2.1).

Nonostante l’ultimo risultato, in generale, non si puo essere sicuri chescegliendo opportuni punti di interpolazione ed un polinomio di interpo-lazione di grado alto, si otterra una buona approssimazione di una certa

1xi = cos( (2i+1)π2(n+1) ), i = 0, . . . , n sono gli n + 1 zeri reali del polinomio di Chebishev

di grado n + 1 definito in [−1, 1]; se i punti di interpolazione sono in [a, b], si applichiun mapping degli zeri da [−1, 1] in [a, b].


n + 1 punti maxx |R(x)| n + 1 punti maxx |R(x)|11 1.92 12 0.5521 58.59 22 17.2931 2277.70 32 665.64

Tabella 2.2: Errore Massimo per l’interpolazione polinomiale dellafunzione di Runge su punti equidistanti.

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Figura 2.1: Interpolazione polinomiale della funzione di Runge su puntiequidistanti; a sinistra 11 punti, a destra 12.

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Figura 2.2: Interpolazione polinomiale della funzione di Runge su puntizeri di Chebishev; a sinistra 11 punti, a destra 12.

funzione. Si puo anche vedere che i polinomi non sono sufficientementeflessibili ad approssimare funzioni che hanno differenti andamenti in diffe-renti parte dell’intervallo di definizione. Questo e il caso di molte funzioni

2.6. INTERPOLAZIONE POLINOMIALE A TRATTI 35

n + 1 punti maxx |R(x)| n + 1 punti maxx |R(x)|11 0.1089 12 0.182821 0.0153 22 0.025331 0.0021 32 0.0035

Tabella 2.3: Errore Massimo per l’interpolazione polinomiale dellafunzione di Runge su punti zeri di Chebyshev.

che descrivono un fenomeno fisico o la forma di un oggetto. Oltre a que-sto si sottolinea il pericolo, dal punto di vista numerico, di lavorare conpolinomi di grado alto. Numericamente parlando si dovrebbero usare solopolinomi di grado basso. Queste considerazioni sono alla base dell’interpo-lazione con polinomi a tratti che consiste nell’approssimare una funzionecon differenti polinomi di grado basso in differenti parti dell’intervallo diinteresse.

2.6 Interpolazione polinomiale a tratti

Con interpolazione polinomiale a tratti si intende l’interpolazione di unset di dati su un intervallo con piu polinomi ciascuno dei quali definitoin un sottointervallo dell’intervallo dato. In particolare siano assegnatem+1 osservazioni in corispondenza di m+1 punti distinti e ordinati, cioesiano assegnati (xi, yi) i = 0, . . . ,m con xi < xi+1. Allora un intepolantepoliniomiale a tratti consiste in m polinomi pi(x), i = 0, . . . ,m − 1 digrado n, definiti sugli intervalli [xi, xi+1] con le seguenti proprieta:

•pi−1(xi) ≡ pi(xi) = yi i = 1, . . . ,m− 1

• fissato k ≤ n− 1, non negativo, deve valere

p(ℓ)i−1(xi) ≡ p

(ℓ)i (xi) ℓ = 1, . . . , k i = 1, . . . ,m− 1

Si dira che l’interpolante a tratti e di classe Ck intendendo che e unafunzione continua fino alla derivata di ordine k; quando k = n − 1 siha la massima continuita e l’interpolante polinomiale a tratti viene dettafunzione spline.

Prima di procedere a presentare due metodi di interpolazione moltoutilizzati in pratica si richiama l’interpolazione di Hermite con un po-linomio cubico nella base di Bernstein; questa consiste nel determina-re il polinomio cubico p(t) che soddisfa le seguenti quattro condizioni di


interpolazione:

p(0) = y0, p′(0) = y′0, p(1) = y1, p′(1) = y′

1.

Si dovranno determinare i coefficienti b0, b1, b2 e b3. Due di loro sonobanalmente determinati:

b0 := y0, b3 := y1.

Per i rimanenti, richiamiamo il valore della derivata di un polinomio nellabase di Berstein negli estremi (osservazione 1.5):

p′(0) = 3(b1 − b0), p′(1) = 3(b3 − b2).

Da queste si possono ricavare b1 e b2:

b1 := y0 +1

3y′

0, b2 := y1 −1

3y′

1.

2.6.1 Interpolazione locale

Per interpolazione locale si intende un metodo che determina il polinomioa tratti interpolante ricavando ogni singolo polinomio pi(x) indipendente-mente dagli altri. Vediamo questo modo di procedere nel caso particolaredi polinomi a tratti cubici con regolarita C1; i singoli polinomi cubicipi(x) possono essere espressi nella base di Bernstein e determinati in mo-do da interpolare le due osservazioni yi e yi+1 in corrispondenza di xi exi+1 sull’intervallo [xi, xi+1] e i valori di derivata y′

i e y′i+1 sempre in cor-

rispondenza di xi e xi+1. Si tratta di interpolanti cubici di Hermite chebanalmente costituiranno un unico interpolante cubico a tratti dei punti(xi, yi) i = 0, . . . ,m e di classe C1; la Fig. 2.3 mostra il risultato di unatale interpolazione sulla funzione test di Runge su punti equidistanti.

N. punti maxx |R(x)| N. punti maxx |R(x)|11 0.0182 12 0.111421 0.0111 22 0.018131 0.0042 32 0.0048

Tabella 2.4: Errore Massimo per l’interpolazione polinomiale cubica atratti C1 della funzione di Runge su punti equidistanti e derivate stimate.

Operativamente ogni intervallo [xi, xi+1] viene idealmente mappato nel-l’intervallo [0, 1] ed in questo si applica una interpolazione cubica di Hermi-te vista, nella base di Bernstein. Le derivate assegnate y′

i dovranno essere


−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Figura 2.3: Interpolazione polinimiale cubica a tratti C1 della funzionedi Runge su punti equidistanti e derivate stimate; a sinistra 11 punti, adestra 12.

N. punti maxx |R(x)| N. punti maxx |R(x)|11 0.0129 12 0.029321 0.0013 22 0.002931 0.0005 32 0.0006

Tabella 2.5: Errore Massimo per l’interpolazione polinomiale cubicaa tratti C1 della funzione di Runge su punti equidistanti e derivateanalitiche.

opportunamente scalate e per la precisione le osservazioni da interpolaresaranno:

(0, yi) (0, y′ihi+1) (1, yi+1) (1, y′

i+1hi+1)

con hi = xi+1 − xi i = 0, . . . ,m− 1.

Stima delle derivate

I metodi di interpolazione locale fanno uso delle derivate y′i in corrispon-

denza dei punti assegnati xi. Se queste non sono date, devono esserecalcolate come parte del problema di interpolazione. In letteratura esisto-no un certo numero di metodi per ottenere stime delle derivate a partiredai dati assegnati. Qui ne riportiamo una semplice classe: siano

hi = xi+1 − xi i = 0, . . . ,m− 1fi = yi+1 − yi i = 0, . . . ,m− 1

di = fi

hi

i = 0, . . . ,m− 1


allora si definisce una stima delle derivate come:

Di = (1− αi)di−1 + αidi i = 1, . . . m− 1.

A seconda di come vengono scelti i parametri αi si hanno differenti schemi;la scelta

αi =hi−1

hi−1 + hi

i = 1, . . . m− 1

porta alla stima di Bessel. Si noti che nel primo ed ultimo punto le sti-me date non sono definite e che questi devono essere trattati come casiparticolari. Per la stima di Bessel si definisce:

D0 = 2d0 −D1 Dm = 2dm−1 −Dm−1.

2.6.2 Interpolazione globale

Per interpolazione globale si intende un metodo che determina il polinomioa tratti interpolante ricavando i polinomi pi(x) in modo dipendente daglialtri. Vediamo questo modo di procedere nel caso particolare di polinomia tratti cubici con massima regolarita C2, anche detta spline cubica diinterpolazione.Esprimiamo i singoli polinomi come polinomi di Bernstein in [xi, xi+1] concoefficienti qi, ri, si e qi+1 per i = 0, . . . ,m − 1. Affinche due polino-mi consecutivi, nella base di Bernstein, si raccordino C2 devono esseresoddisfatte le seguenti condizioni o insieme di condizioni:

qi − si−1

hi−1

=ri − qi

hi

i = 1, . . . ,m− 1

qi − 2si−1 + ri−1

h2i−1

=si − 2ri + qi

h2i

i = 1, . . . ,m− 1.(2.2)

Si tratta di 2m − 2 equazioni lineari nelle 2m incognite ri ed si, i =0, . . . ,m − 1 (i qi, per soddisfare le condizioni di interpolazione, devonovalere yi, i = 0, . . . ,m); solitamente si aggiungono due ulteriori condizionio si assegnano liberamente i valori per due incognite (frequentemente perr0 ed sm−1), cosı da avere un sistema quadrato. Il dover risolvere tale siste-ma lineare equivale a determinare i polinomi a tratti in modo dipendentegli uni dagli altri.

Osservazione 2.2 Le condizioni sopra date derivano dalla definizione dipi(x) i = 0, . . . ,m − 1 e dalle espressioni delle loro derivate prima eseconda;

pi(x) = qiB0,3(x) + riB1,3(x) + siB2,3(x) + qi+1B3,3(x)


p′i(x) =3

hi

[(ri − qi)B0,2(x) + (si − ri)B1,2(x) + (qi+1 − si)B2,2(x)]

p′′i (x) =6

h2i

[(si − 2ri + qi)B0,1(x) + (qi+1 − 2si + ri)B1,1(x)]

Se poniamo

Di =qi − si−1

hi−1

=ri − qi

hi

(2.3)

possiamo riscrivere le seconde condizioni 2.2 solo in termini di Di e qi so-stituendo al posto di si ed ri le espressioni trovate mediante la 2.3; sara:

qi−1 − qi + 2Dihi−1 + Di−1hi−1

h2i−1

=qi+1 − qi −Di+1hi − 2Dihi

h2i

;

moltiplicando entrambi i membri per hihi−1 e raccogliendo opportunamen-te i termini si ottiene:

hiDi−1 + 2Di(hi + hi−1) + hi−1Di+1 =hi

hi−1

(qi − qi−1) +hi−1

hi

(qi+1 − qi)

per i = 1, . . . ,m− 1.Si tratta di un sistema di m − 1 equazioni nelle m + 1 incognite Di; sipossono allora assegnare D0 e Dm (derivate agli estremi) e risolvere ilseguente sistema quadrato (m − 1) × (m − 1) che, essendo la matrice adiagonale dominante, avra una ed una sola soluzione

2(h1 + h0) h0 0 . . . 0h2 2(h2 + h1) h1 . . . 0. . .0 . . . 0 hm−1 2(hm−1 + hm−2)

D1

D2...Dm−1

=

=

h1

h0

(q1 − q0) + h0

h1

(q2 − q1)− h1D0h2

h1

(q2 − q1) + h1

h2

(q3 − q2)...

hm−1

hm−2

(qm−1 − qm−2) + hm−2

hm−1

(qm − qm−1)− hm−2Dm

.

Determinate le D1, D2, . . . Dm−1 insieme alle assegnate D0 e Dm si puoprocedere al calcolo degli m polinomi che costituiscono la polinomialecubica a tratti di interpolazione. Tali polinomi possono essere valutati in[0, 1] facendo uso dei coefficienti:

qi = yi i = 0, . . . ,m

qi, qi + hiDi, qi+1 − hiDi+1, qi+1 i = 0, . . . ,m− 1.


Osservazione 2.3 Si noti che nel caso i punti siano tali per cui hi ≡1 ∀i, il sistema si semplifica in

4 1 0 0 . . . 01 4 1 0 . . . 00 1 4 1 . . . 0. . .0 . . . 0 0 1 4

D1

D2

D3...Dm−1

=

q2 − q0 −D0

q3 − q1

q4 − q2...

qm − qm−2 −Dm

.

La Fig. 2.4 mostra il risultato dell’interpolazione della funzione test diRunge con una spline cubica su punti equidistanti.

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Figura 2.4: Interpolazione polinomiale cubica a tratti C2 della funzione diRunge su punti equidistanti; a sinistra 11 punti, a destra 12.

N. punti maxx |R(x)| N. punti maxx |R(x)|11 0.02193 12 0.0838221 0.00318 22 0.0080231 0.00084 32 0.0013141 0.00063 42 0.00061

Tabella 2.6: Errore Massimo per l’interpolazione polinomiale cubica atratti C2 della funzione di Runge su punti equidistanti.

Un esempio: controllo di un robot

Le curve spline sono utilizzate in numerose applicazioni industriali e sonoalla base dei moderni sistemi CAD. In questo paragrafo vogliamo vedere


un semplice esempio applicativo di una spline cubica di interpolazione.Nella robotica e frequente la situazione di avere dei motori che muovonodei bracci meccanici i cui movimenti vengono gestiti da un calcolatoreper compiere delle azioni ripetitive. Per focalizzare si puo pensare adun semplice braccio meccanico, come quello schematizzato in Fig. 2.5,composto da due parti (la spalla di lunghezza L1 e il gomito di lunghezzaL2) rispettivamente mosse da due motori. Quando questi due motoriruotano, il braccio puo posizionare la mano (la parte estrema del gomito)in un punto desiderato del piano su cui il braccio e posto. Il robot e gestitoda un calcolatore che riceve l’input dell’utente (per esempio, la posizionedesiderata della mano). Il computer calcola la rotazione che ogni motoredeve effettuare per spostare la mano dalla posizione iniziale a quella finale.Successivamente, a intervalli regolari di tempo, il computer invia i comandidi rotazione ai due motori. In applicazioni come queste il computer usa unainterpolazione spline per generare i comandi di rotazione. Il vantaggio diuna spline e di fornire un movimento regolare ai motori del robot, evitandobruschi cambi di velocita e/o accelerazione, se regolare fino alla derivataseconda (nel caso cubico la spline e C2).

Possiamo, grazie alla trigonometria, definire le espressione per le coor-dinate del gomito e della mano:

gomito:

xgomito = L1 cos(θ1)ygomito = L1 sin(θ1)

θ

θL

L

1

1

2

2

Figura 2.5: Schema di un braccio meccanico con due motori


mano:

xmano = xgomito + L2 cos(θ1 + θ2) = L1 cos(θ1) + L2 cos(θ1 + θ2)ymano = ygomito + L2 sin(θ1 + θ2) = L2 sin(θ1) + L2 sin(θ1 + θ2)

Queste espresioni permettono di impostare gli angoli del gomito e dellaspalla che sono necessari per posizionare la mano nel punto specificatodalle coordinate (xmano, ymano). Per la precisione si ha:

R2 = x2mano + y2

mano

cos(θ2) =R2 − L2

1 − L22

2L1L2

cos(β) =R2 − L2

1 − L22

2L1R

α = arctan(ymano

xmano

)

θ1 =

{

α + β se θ2 < 0α− β se θ2 ≥ 0

Se si utilizza un polinomio di grado 3 per esprimere gli angoli dei giunti infunzione del tempo per 0 ≤ t ≤ T , e posibile specificare quattro condizioniper ciascun angolo. Due di queste condizioni richiedono che il polinomiopassi per il valore θ(0) e per il valore θ(T ). Le altre due condizioni ri-chiedono che la pendenza del polinomio debba essere nulla all’inizio e allafine, in quanto il robot deve partire da una posizione di riposo e finire inuna posizione di riposo.

Esempio 2.3 Supponiamo che il braccio del robot abbia lunghezza L1 =4m ed L2 = 3m. Supponiamo inoltre che la mano debba spostarsi in 2secondi dal punto x(0) = 6.5, y(0) = 0.0 al punto x(2) = 0.0, y(2) = 6.5;dalle espressioni otteniamo:

θ1(0) = −18.717o θ1(2) = 71.283o

θ2(0) = θ2(2) = 44.049o

I polinomi di interpolazione questi valori e derivate nulle in 0 e 2 sono:

θ1(t) = −18.717 + 67.5t2 − 22.5t3

θ2(t) = 44.049

Se il computer utilizza questi polinomi per generare i comandi da inviareai motori, si nota che la mano compie una traiettoria circolare in quantol’angolo θ2 fra i due bracci resta costante per l’intero percorso, mentrevaria solo θ1.


Osservazione 2.4 Si noti che la traiettoria della mano non e una li-nea retta che congiunge la posizione iniziale e quella finale e che questomovimento potrebbe provocare delle collisioni con gli oggetti vicini.

In molte aplicazioni e necessario controllare la traiettoria della mano conmaggior precisione. Supponiamo di voler ottenere una traiettoria oppor-tuna che vogliamo specificare con una serie di posizioni intermedie chedevono essere raggiunte dalla mano del robot, in istanti unitari di tempo,partendo dalla posizione iniziale e terminando in quella finale.Per ottenere un movimento regolare, non tanto della mano o traiettoria,

−8 −6 −4 −2 0 2 4 6 8

−8

−6

−4

−2

0

2

4

6

8

Traiettoria della mano del robot

x (m)

y (

m)

Figura 2.6: Schema di un braccio meccanico con due motori: output diun programma esempio

quanto dei motori, per evitare bruschi cambi di velocita e accelerazione,si devono utilizzare delle funzioni spline per interpolare i valori di θ1(ti) eθ2(ti) in corrispondenza degli istanti ti prefissati; θ1(ti) e θ2(ti) corispon-dono agli angoli che permettono di posizionare la mano nelle posizioni(x(ti), y(ti)) assegnate agli istanti ti. La Fig.2.6 mostra l’ultimo fotogram-mo generato da un programma di esempio che simula il comportamentodi una tale robot che nel suo movimento con accelerazione continua de-ve partire da una posizione iniziale ed arrivare ad una finale passando in


istanti unitari di tempo attraverso quattro posizioni predefinite. La coro-na circolare rappresentata in figura indica il raggio di azione possibile delrobot avendo fissato per tale simulazione L1 = 6 ed L2 = 3.

2.7 Un’applicazione: curve cubiche a tratti

di interpolazione

Nei pacchetti grafici di disegno 2D, nei sistemi CAD e in numerose altreapplicazioni la determinazione di curve mediante interpolazione di puntirisulta una tecnica molto usuale. In questa sezione si vogliono presentarei due tipi sicuramente piu frequentemente usati; l’interpolazione con unacurva cubica a tratti C1 e con una spline cubica.Siano dati m + 1 punti Qi ≡ (xi, yi)i=0,...,m; si vogliono determinare lacurva cubica a tratti C1 e la spline cubica che interpolano i punti Qi; sicerca una curva in forma parametrica

C(t) =

(

C1(t)C2(t)

)

tale che siano verificate le condizioni di interpolazione:

{

C1(ti) = xi

C2(ti) = yii = 0, . . . ,m

Il problema di interpolazione con una curva (funzione vettoriale a duecomponenti) si riduce a due problemi di interpolazione con funzioni scalari,cioe si cercano: la funzione C1(t) che interpola i punti (ti, xi) e la funzioneC2(t) che interpola i punti (ti, yi).Si noti che i valori parametrici ti in corrispondenza dei quali interpolarenon sono assegnati e devono essere determinati come parte del problemastesso. Vediamo due possibili modi.

Metodo della parametrizzazione uniforme. Consiste nel suddivide-re l’intervallo [0, 1] in m sottointervalli di uguali ampiezze, ovveronel prendere m + 1 punti equidistanti:

ti =i− 1

mi = 1, . . . ,m.

Questo tipo di parametrizzazione, pero, non tiene conto della distan-za relativa dei punti sulla curva.Per ovviare a questo inconveniente si usa il metodo seguente:

2.7. UN’APPLICAZIONE: CURVE DI INTERPOLAZIONE 45

Metodo della parametrizzazione della corda. In questo caso si tie-ne conto della geometria dei punti, poiche si mantengono invariati irapporti:

ti+1 − titi+2 − ti+1

=‖Qi+1 −Qi‖2‖Qi+2 −Qi+1‖2

.

Per determinare i valori ti si puo procedere come segue:

τ1 = 0

τi = τi−1 +√

(xi − xi−1)2 + (yi − yi−1)2 i = 2, . . . ,m + 1

Infine si ottengono i valori:

ti =τi

τn

i = 1, . . . ,m + 1.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figura 2.7: Interpolazione con curva polinomiale cubica a tratti C1 constima dei vettori derivati nei punti (sinistra); vettori derivati scalati di 1/2(centro); vettori derivati scalati di 1/4 (destra).

Determinati i parametri ti, i metodi piu utilizzati sono quelli visti nellesezioni precedenti e per la precisione l’interpolazione polinomiale cubicaa tratti C1 e spline cubica che fornisce una curva C2. La prima e moltointeressante per il bassissimo costo computazionale e per la possibilita diavere m+1 parametri di forma consistenti nei moduli dei vettori tangentinei punti di interpolazione; tali scalari giocano come parametri di tensio-ne locale per la curva o globali se modificati tutti contemporaneamente.


−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

Figura 2.8: Interpolazione con curva spline cubica di 33 punti.

La Fig.2.7 mostra tre interpolazioni di 8 punti disposti a forma di essecon cubiche a tratti C1 dove i vettori derivati, inizialmente stimati comeindicato nella sezione 2.6.1, sono stati globalmente scalati di 1/2 e poi1/4 evidenziando il tensionamento relativo. La Fig.2.8 mostra invece unacurva di interpolazione spline cubica con parametrizzazione della corda di32 punti disposti nel piano a definire un profilo 2D di Paperino.

Capitolo 3

Approssimazione ai minimiquadrati

Siano assegnate m osservazioni in corispondenza di m punti distinti, cioesiano assegnati (xi, yi) i = 1, . . . ,m. Nel caso generale viene individuatouno spazio Φ a dimensione finita n + 1 ≤ m, di funzioni reali e una lorobase di rappresentazione φ0, φ1, . . . , φn cosı che ogni elemento φ ∈ Φ siadefinita da una loro combinazione lineare, cioe

φ =n∑

i=0

aiφi(x).

Il problema della migliore approssimazione nel senso dei minimi quadra-ti consiste nel determinare una funzione funzione φ∗ ∈ Φ (od anche icoefficienti a∗

0, a∗1, . . . , a

∗n che la definiscono) tale che

m∑

k=1

[φ∗(xk)− yk]2 ≤

m∑

k=1

[φ(xk)− yk]2 ∀φ ∈ Φ.

Questo equivale a determinare il minimo della seguente funzione in n + 1variabili:

g(a0, a1, . . . , an) =m∑

k=1

[n∑

i=0

aiφi(xk)− yk]2

cioeminai∈R

g(a0, a1, . . . , an).

Un metodo per risolvere questo problema di minimo consite nell’impostaree risolvere il sistema delle equazioni normali, che deriva dall’imporre

∂g(a0, a1, . . . , an)

∂ai

= 0 i = 0, . . . , n

47

48 CAPITOLO 3. APPROSSIMAZIONE AI MINIMI QUADRATI

o in forma esplicita

∂

∂ai

m∑

k=1

[

n∑

i=0

aiφi(xk)− yk

]2

= 0 i = 0, . . . , n.

Infatti, essendo la funzione g quadratica con coefficienti positivi per itermini a2

i , avra un minimo. Pima di procedere alla derivazione riscriviamola funzione g sviluppando il quadrato:

g(a0, a1, . . . , an) =m∑

k=1

y2k−2

n∑

i=0

ai

m∑

k=1

ykφi(xk)+n∑

i=0

n∑

j=0

aiaj

m∑

k=1

φi(xk)φj(xk)

allora

∂g

∂ai

= −m∑

k=1

ykφi(xk) +n∑

j=0

aj

(

m∑

k=1

φi(xk)φj(xk)

)

i = 0, . . . , n

da cui il sistema lineare

Ga = r

dove

G = {gi,j}i,j=0,...,n con gi,j =m∑

k=1

φi(xk)φj(xk)

ed

r = {ri}i=0,...,n con ri =m∑

k=1

ykφi(xk).

Osservazione 3.1 Si noti che algoritmicamente, il sistema Ga = r puoessere costruito calcolando la matrice H = {hi,j}i=1,...,m,j=0,...,n delle fun-zioni base nei punti, cioe hi,j = φj(xi) e quindi ottenere la matrice G edil vettore r come:

G := HT H r := HTy.

Nel caso generale qui descritto, la scelta dello spazio Φ e solitamente det-tato dal fatto che esista sempre una soluzione unica al problema, dallaregolarita desiderata per l’approssimante e dalla particolarita dei dati.

In queste dispense si trattera solamente il caso dell’approssimazionepolinomiale, in cui la matrice G risulta simmetrica e definita positiva, ilche garantisce sempre l’esistenza e l’unicita della soluzione ed una buonaregolarita della stessa.

3.1. FORMA MONOMIALE 49

3.1 Forma Monomiale

Teorema 3.1 Siano dati m punti (xi, yi), i = 1, . . . ,m con xi distinti, al-lora esiste ed e unico il polinomio p ∈ Pn, con n ≤ m, di approssimazionenel senso dei minimi quadrati.

dim. Si consideri p(x) nella forma monomiale e si noti la forma dellamatrice H:

H =

1 x1 x21 . . . xn

1

1 x2 x22 . . . xn

2

1 x3 x23 . . . xn

3

. . .

. . .1 xm x2

m . . . xnm

si tratta di una matrice di Vandermonde rettangolare, ed essendo i puntixi distinti, sara di rango massimo n + 1; segue che la matrice G = HT Hsara non singolare e che il sistema lineare delle equazioni normali avraun’unica soluzione; segue che p(x) di migliore approssimazione esiste ed eunico.

3.2 Forma di Bernstein

Nella sezione precedente si e dimostrata l’esistenza ed unicita del polino-mio approssimante nella base monomiale, ma tale polinomio puo essererappresentato in una qualunque altra base. Procediamo quindi alla suadeterminazione via soluzione di un sistema lineare, ma a partire dalla basedi Bernstein. Siano dati m punti (xi, yi), i = 1, . . . ,m con xi distinti in[a, b] e si vuole determinare p(x) di grado n ≤ m, nella base di Bernsteintale che risulti il polinomio di approssimazione nel senso dei minimi qua-drati. Siano (ti, yi), i = 1, . . . ,m i punti traslati e scalati in [0, 1] mediantela 1.6 e sia

p(t) =n∑

i=0

biBi,n(t) con t ∈ [0, 1]

la forma polinomiale che vogliamo approssimi i valori assegnati. Costruen-do la marice H delle funzioni di Bernstein nei punti si puo impostare erisolvere il sistema delle equazioni normali. Il vettore soluzione b delsistema lineare fornira sia il polinomio approssimante in [0, 1] che in [a, b].

Si noti che la matrice G delle equazioni normali, nel caso si usi la basedi Bernstein nei punti, risulta meglio condizionata che non la matrice Gnel caso si usi la forma monomiale (si veda Tab.3.1).


grado n Bernst. Potenze10 2.77E+05 1.01E+1415 2.73E+08 6.12E+2120 3.95E+11 1.96E+2825 9.68E+15 1.95E+3630 2.18E+16 3.64E+3935 1.75E+17 1.15E+51

Tabella 3.1: Indici di condizionamento delle matrici G ottenute valutandoi polinomi di Bernstein e la base delle potenze in 51 punti equidistanti.

3.3 Forma ortogonale

Come gia detto sia nell’introduzione che nel capitolo dell’interpolazione,a seconda del problema, si puo determinare una base piu idonea delle al-tre sia dal punto di vista del condizionamento, ma anche che permetta diusare un metodo piu efficiente oltre che stabile. Nel caso dell’approssi-mazione polinomiale sicuramente la base piu idonea e quella dei polinomi{pi(x)} i = 0, . . . , n ortogonali sui punti dati {xk} k = 1, . . . ,m, ossiacaratterizzati dalla proprieta:

m∑

k=1

pi(xk)pj(xk) = 0 i 6= j.

Con tale base la matrice G delle equazioni normali risulta diagonale e lasoluzione e semplicemente data da:

ai =

∑mk=1 ykpi(xk)

∑mk=1[pi(xk)]2

i = 0, . . . , n

Per costruire tale base ortogonale di polinomi di grado n su m punti {xk}assegnati, si puo usare lo schema ricorrente proposto da Hoseholder eStiefel.

3.3.1 Generazione polinomi ortogonali

Vediamo i polinomi ortogonali su un insieme discreto di punti, suggeritida Householder e Stiefel.

{

p0(x) = 1pi(x) = xpi−1(x)− αipi−1(x)− βipi−2(x) i = 1, 2, . . .

3.4. ESEMPI NUMERICI 51

dove

αi =

∑mk=1 xk[pi−1(xk)]

2

wi−1

i = 1, 2, . . .

β1 = 0

βi =wi−1

wi−2

i = 2, 3, . . .

con

wi =m∑

k=1

[pi(xk)]2 i = 0, 1, . . . .

EsempioSiano assegnati in [0, 1] i cinque punti xi = (i − 1)h i = 1, . . . , 5 conh = 1/4. Allora i polinomi ortogonali su tali punti che generano P1 sono:

{

p0(x) = 1p1(x) = xp0(x)− α1p0(x)

dove

α1 =

∑5k=1 xk[p0(xk)]

2

w0

=0 + 0.25 + 0.5 + 0.75 + 1

5= 0.5

essendo

w0 =5∑

k=1

[p0(xk)]2 = 5

da cuip1(x) = x− 0.5

Verifichiamo che5∑

k=1

p0(xk)p1(xk) = 0;

infatti

1 · (−0.5) + 1 · (−0.25) + 1 · 0 + 1 · 0.25 + 1 · 0.5 = 0.

3.4 Esempi numerici

Le Fig.3.1 e 3.2 mostrano il risultato di un’approssimazione ai minimi qua-drati di 51 punti equidistanti della funzione di Runge rispettivamente conpolinomi di grado 10, 15, 30 e 35 nella base monomiale. I differenti indi-ci di condizionamento della Tab.3.1 indicano una elaborazione numericapotenzialmente meno accurata nel caso monomiale che non nella base di


−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Figura 3.1: Approssimazione polinomiale nella base delle potenze dellafunzione di Runge di 51 punti equidistanti; a sinistra grado 10, a destra15.

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Figura 3.2: Approssimazione polinomiale nella base delle potenze dellafunzione di Runge di 51 punti equidistanti; a sinistra grado 30, a destra35.

bernstein; in realta questa differenza non risulta graficamente apprezza-bile ne usando Bernstein ne usando i polinomi ortogonali. Determinatoil polinomio p∗(x) di grado n di migliore approssimazione, chiameremoresiduo la quantita

dn =n∑

i=0

[p∗(xk)− yk]2.

Si nota che all’aumentare di n (n < m), sara dn+1 < dn, e per n = m− 1,dm−1 ≡ 0: questo, come e ben visibile dalle Fig.3.2, non assicura che p∗(x)converga alla funzione test ne analiticamente ne numericamente. Anche in

3.5. RETTA DI REGRESSIONE LINEARE 53

questo caso i polinomi risultano poco flessibili e numericamente instabiliper gradi alti.

Per l’approssimazione nel senso dei minimi quadrati e possibile utiliz-zare polinomi a tratti o spline tipicamente di grado basso, con il vantaggiodi avere delle ottime funzioni di approssimazione senza i problemi numericie di flessibilita tipici dei polinomi di grado alto.

3.5 Retta di regressione lineare

Nelle applicazioni pratiche, noto un set di dati sperimentali acquisiti daun fenomeno di andamento lineare, e molto utilizzata la tecnica nota comeretta di regressione lineare che consiste nel determinare il polinomiolineare di approsimazione nel senso dei minimi quadrati dei dati assegnati.La determinazione di tale polinomio viene effettuata mediante la soluzionedel sistema lineare delle equazioni normali, che in questo caso risultera didimensioni 2 × 2 e la cui soluzione puo essere data in forma esplicitarisolvendo analiticamente il sistema suddetto. Sara:

G =

(

∑mk=1 1

∑mk=1 xk

∑mk=1 xk

∑mk=1 x2

k

)

e definendo

x =1

m

m∑

k=1

xk y =1

m

m∑

k=1

yk

il baricentro dei punti dati, il sistema puo essere scritto in modo semplifi-cato come:

(

m mxmx

∑mk=1 x2

k

)(

a0

a1

)

=

(

my∑m

k=1 xkyk

)

da cui con semplici passaggi si ricava

a0 = y − a1x a1 =

∑mk=1 xkyk −mxy∑m

k=1 x2k −mx2

(3.1)

In modo ancor piu operativo, l’espressione del polinomio di approssima-zione puo essere scritta come

p(x) = a0 + a1x = y + a1(x− x)

che comporta solo la valutazione di a1 dalla seconda delle 3.1 e geome-tricamente dice che la retta di regressione lineare passa per il baricentro


del set di dati ed ha pendenza a1. In questo caso diremo che la rettadi regressione lineare e rappresentata nella base con centro x data dallefunzioni {1, (x− x)}.

Se richiesto, la retta di regressione lineare si puo esprimere semplice-mente nella base lineare di Bernstein in [0, 1] come

p(t) = a0(B0,1(t) + B1,1(t)) + a1B1,1(t) =

= a0B0,1(t) + (a0 + a1)B1,1(t).

3.6 Un’applicazione: curve di approssima-

zione

Possiamo risolvere problemi di approssimazione ai minimi quadrati confunzioni vettoriali o curve polinomiali o polonomiali a tratti, per esempionel piano.Siano assegnati m punti nel piano euclideo Qi = (xi, yi)i=1,...,m. Il nostroobiettivo e determinare una curva polinomiale nella base di Bersntein:

C(t) =

(

C1(t)C2(t)

)

=n∑

i=0

PiBi,n(t)

che rende minima la seguente espressione:

m∑

i=1

‖Qi − C(ti)‖22 =

m∑

i=1

((xi − C1(ti))2 + (yi − C2(ti))

2)

Come nel caso dell’interpolazione, il problema della determinazione diuna curva polinomiale che approssimi i punti Qi, puo essere scompostoin due sottoproblemi, cioe nella ricerca di due funzioni polinomiali C1(t)e C2(t) che approssimino nel senso dei minimi quadrati rispettivamente ipunti (ti, xi) e (ti, yi) per i = 1, . . . ,m. Risulta sufficiente trovare le duefunzioni componenti della curva applicando due volte il metodo dei minimiquadrati.Il primo passo da fare e, ancora una volta, la determinazione dei parametriti a cui far corrispondere i punti Qi. Metodi utilizzabili sono, come giavisto nel caso dell’interpolazione con curve, la parametrizzazione uniformee la parametrizzazione della corda.

Capitolo 4

Approssimazione di forma

In questa sezione si vuole mostrare come i polinomi di Bernstein, inizial-mente introdotti da Bernstein stesso per l’approssimazione uniforme su unintervallo, godono anche della proprieta di essere approssimanti di formao in modo piu tecnico di essere approssimanti variation diminishing (VD).

4.1 Approssimazione uniforme

I polinomi di Bernstein furono originariamente introdotti nell’approssima-zione di funzioni continue f(x) su un intervallo [a, b]. Dati n + 1 puntidi f(x) equidistanti in [a, b], l’approssimazione di Bernstein di grado n,Bn[f(x)] di f(x) e definita come

Bn[f(x)] =n∑

i=0

f(ξi)Bi,n(x)

con ξi = f(a + i(b − a)/n) e Bi,n(x) i polinomi di Bernstein di gradon in [a, b]. Questa approssimazione puo essere costruita soddisfando unqualunque errore δ usando polinomi di grado sufficientemente elevato, cioeesiste un valore nδ tale che:

|Bn[f(x)]− f(x)| < δ

∀x ∈ [a, b] quando n > nδ. Cosı si dice che l’approssimazione di Bernsteinconverge uniformemente ad f(x):

limn→∞

Bn[f(x)] = f(x) x ∈ [a, b].

Questa proprieta della base di Bernstein permette di dimostrare il teore-ma di Weierstrass che una funzione continua su un intervallo puo essereapprossimata con una data tolleranza da un polinomio di opportuno gradon.

55

56 CAPITOLO 4. APPROSSIMAZIONE DI FORMA

4.2 Approssimazione VD

Definizione 4.1 Variazione di segno di un vettoreSia c = (c0, c1, . . . , cn) un vettore ad elementi reali; si definisce numero divariazioni di segno di c in senso forte, e lo si indica con

V −[c]

il numero di variazioni di segno nella sequenza c0, c1, . . . , cn dove si igno-rano gli zeri.

Esempio 4.1 Sia c = (−1, 3.2, 1.5,−3.7, 4) allora V −[c] = 3.

Definizione 4.2 Variazione di segno di una funzioneSia f una funzione reale in [a, b]; si dira che:

V −[f(z)] = supn+1{V−[f(z0), . . . , f(zn)] con z0 < z1 < . . . < zn}.

Sia data una funzione f(x) almeno continua nell’intervallo [a, b]. Nelcaso generale viene individuato uno spazio Φ a dimensione finita n + 1, difunzioni reali e una loro base di rappresentazione φ0, φ1, . . . , φn cosı cheogni elemento φ ∈ Φ e definito da una loro combinazione lineare, cioe

φ =n∑

i=0

aiφi(x).

Il problema dell’approssimazione di forma o variation diminishing con-siste nel determinare una funzione φ∗ ∈ Φ tale che

1. se f(x) e lineare, allora la φ∗ e lineare e coincide con la f(x);

2.V −[φ∗(x)] ≤ V −[f(x)] x ∈ [a, b]

che significa che la variazione in segno in senso forte della φ∗ e minoreo uguale alla variazione in segno in senso forte della f(x).

Se queste due proprieta vengono soddisfatte si ha come conseguenza che

V −[φ∗(x)− (a0 + a1x)] ≤ V −[f(x)− (a0 + a1x)] x ∈ [a, b] (4.1)

che dice che la φ∗(x) non puo essere piu oscillante della f(x) che approssi-ma, in quanto il numero di intersezioni della φ∗(x) con una retta e sempreminore o uguale al numero di intersezioni della f(x) con quella retta.

4.2. APPROSSIMAZIONE VD 57

Teorema 4.1 Sia n ≥ 1, allora ∀x ∈ [a, b] si ha:

n∑

i=0

ξiBi,n(x) = x

con ξi = a+i(b−a)n

.

Teorema 4.2 Sia p : [a, b]→ R un polinomio nella base di Bernstein e siconsideri come approssimazione di una f(x) : [a, b]→ R almeno continua,la Bn[f(x)] definita nella sezione 4.1; allora Bn[f(x)] e approssimante VDdella f(x).

dim.Verifichiamo le due proprieta caratterizzanti:1)Bn[f(x)] deve preservare le funzioni lineari:

Bn[a0 + a1x] =n∑

i=0

(a0 + a1ξi)Bi,n(x) =

= a0

n∑

i=0

Bi,n(x) + a1

n∑

i=0

ξiBi,n(x) = a0 + a1x

dove nella prima sommatoria si e applicata la proprieta di partizione del-l’unita mentre nella seconda il Teorema 4.1.2) La Bn[f(x)] gode della proprieta di variazione in segno dei coefficienti(vedi proprieta sui polinomi di Bernstein), cioe:

V −[Bn[f(x)]] ≤ V −[f(ξ0), . . . , f(ξn)]

poiche gli ξi sono in ordine crescente, individueranno sicuramente unan + 1-pla che soddisfa la definizione 4.2 di V −[f(x)] per cui sara

V −[f(ξ0), . . . , f(ξn)] ≤ V −[f(x)]

verificando cosı la seconda proprieta.

La Fig.4.1 mostra una funzione f lineare a tratti e alcuni polinomidi approssimazione VD di grado crescente. Dai grafici si vede come alcrescere del grado i polinomi approssimano sempre meglio la funzione f ;vale infatti il teorema di Weierstrass che la Bn[f(x)] verifica. La stessafigura mostra anche come ogni polinomio approssima in forma, secondo laproprieta conseguente 4.1, la funzione data.


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−0.5

0

0.5

1

1.5

n=3

n=6

n=12

n=24

n=48

Figura 4.1: Approssimazione polinomiale VD di una funzione lineare atratti.

4.3 Significato geometrico dei coefficienti

Dato un polinomio di grado n nella base di Bernstein, i coeffcienti bi, tipi-camente scalari, assumono un importante significato geometrico derivantedall’approssimazione VD vista. Infatti, tali coefficienti possono essere in-terpretati come i valori di una funzione continua in corrispondenza delleascisse ξi i = 0, . . . , n che il polinomio approssima nel senso VD. Per ana-logia con il caso di funzioni vettoriali (curve di Bezier) la funzione continuaviene scelta come la lineare a tratti di vertici i punti (ξi, bi).E noto che una funzione scalare puo essere sempre rappresentata come unafunzione vettoriale in modo banale. Questo si puo applicare ad un polino-mio nella base di Bernstein e rappresentarlo come una funzione vettoriale,ed in particolare per il Teorema 4.1 come una curva di Bezier.

Sia p(x) =∑n

i=0 biBi,n(x) con x ∈ [a, b], un polinomio scalare nella basedi Bernstein. Allora questo puo essere rappresentato in forma vettorialecome:

C(t) =

(

t∑n

i=0 biBi,n(t)

)

e per il teorema 4.1

C(t) =

(

∑ni=0 ξiBi,n(t)

∑ni=0 biBi,n(t)

)

.

4.3. SIGNIFICATO GEOMETRICO DEI COEFFICIENTI 59

Questa rappresentazione giustifica, che nel caso scalare, la funzione chep(x) approssima nel senso VD sia assunta come la lineare a tratti di verticii punti (ξi, bi), i = 0, . . . , n.

Osservazione 4.1 Dalla Fig.4.1 si puo notare che la convergenza del po-linomio approssimante VD e molto lenta, ossia si ottiene per polinomidi grado molto elevato. Anche in questo caso si puo ricorrere all’uso dipolinomi a tratti o spline tipicamente di grado basso (cubiche) e ottenereun’approsimazione VD che converge alla f(x) ∈ C[a,b] all’aumentare deitratti o intervalli in cui si partiziona [a, b].

Elenco delle figure

1.1 Interpolazione di punti 3D con una funzione vettoriale . . . 2

1.2 Amplificazione dell’errore . . . . . . . . . . . . . . . . . . . 8

1.3 Polinomi base di Bernstein di grado 3. . . . . . . . . . . . 10

1.4 Grafici delle valutazioni effettuate in 33 punti. . . . . . . . 14

1.5 Grafici delle valutazioni effettuate in 201 punti. . . . . . . 15

1.6 Esempio di curva di Bezier; caso n = 3 a sinistra; modificadi forma mediante spostamento di un punto di controllo adestra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.1 Interpolazione polinomiale della funzione di Runge su puntiequidistanti; a sinistra 11 punti, a destra 12. . . . . . . . . 34

2.2 Interpolazione polinomiale della funzione di Runge su puntizeri di Chebishev; a sinistra 11 punti, a destra 12. . . . . . 34

2.3 Interpolazione polinimiale cubica a tratti C1 della funzionedi Runge su punti equidistanti e derivate stimate; a sinistra11 punti, a destra 12. . . . . . . . . . . . . . . . . . . . . . 37

2.4 Interpolazione polinomiale cubica a tratti C2 della funzionedi Runge su punti equidistanti; a sinistra 11 punti, a destra12. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

2.5 Schema di un braccio meccanico con due motori . . . . . . 41

2.6 Schema di un braccio meccanico con due motori: output diun programma esempio . . . . . . . . . . . . . . . . . . . . 43

2.7 Interpolazione con curva polinomiale cubica a tratti C1 constima dei vettori derivati nei punti (sinistra); vettori deri-vati scalati di 1/2 (centro); vettori derivati scalati di 1/4(destra). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

2.8 Interpolazione con curva spline cubica di 33 punti. . . . . . 46

3.1 Approssimazione polinomiale nella base delle potenze dellafunzione di Runge di 51 punti equidistanti; a sinistra grado10, a destra 15. . . . . . . . . . . . . . . . . . . . . . . . . 52

61

62 ELENCO DELLE FIGURE

3.2 Approssimazione polinomiale nella base delle potenze dellafunzione di Runge di 51 punti equidistanti; a sinistra grado30, a destra 35. . . . . . . . . . . . . . . . . . . . . . . . . 52

4.1 Approssimazione polinomiale VD di una funzione lineare atratti. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Bibliografia

[BBCM92] R. Bevilacqua, D. Bini, M. Capovani, O. Menchi. MetodiNumerici. Zanichelli, 1992.

[Bez70] P. Bezier. Numerical Control; Mathematics and Applications.Johm Wiley & Sons, 1970.

[FaPr87] I. D. Faux, M. J. Pratt. Computational Geometry for Design andManufacture. John Wiley & Sons, 1987.

[FaRa87] R. T.Farouki, V. T.Ragian. On the numerical condition of po-lynomials in Bernstein form; Computer Aided Geometric Design, 4(1987), 191-216.

[FaRa88] R. T.Farouki, V. T.Ragian. Algorithms for polynomials inBernstein form; Computer Aided Geometric Design, 5 (1988), 1-26.

63

ARGOMENTI DEL CORSO CALCOLO NUMERIC0vvw.web.cs.unibo.it/wiki/images/0/09/CC-2009_2010-02.pdf ·...

Documents

Transcript of ARGOMENTI DEL CORSO CALCOLO NUMERIC0vvw.web.cs.unibo.it/wiki/images/0/09/CC-2009_2010-02.pdf ·...