Chimica Industriale 1

Interpolazione e Approssimazione di Funzioni

Uno dei primi problemi della matematica e anche uno dei piu applicati e la costruzione diuna approssimazione di una funzione data f mediante funzioni piu semplici, tipicamente,ma non sempre, polinomi. Una variante di questo problema e la costruzione di una funzioneregolare a partire da un insieme discreto di punti. Quando si deve trattare una funzionecomplicata a tal punto che non sia possibile trattare i suoi valori in modo diretto si cerca diapprossimarla con una funzione nota e piu semplice da calcolare. Un modo per risolvere ilproblema e di definire una classe di funzioni approssimanti e di cercare la soluzione migliorein funzione di una opportuna metrica. Una scelta possibile e la classe dei polinomi di gradoassegnato m, Pm.Definizione di Polinomio

• Una funzione p definita ∀x ∈ R da:

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

n∑

i=0

aixi

dove n e un intero non negativo e a0, a1, . . . an sono numeri reali fissati, e dettoPolinomio;

• se p(x) ha questa rappresentazione e an 6= 0 allora p(x) ha grado n (ordine n+1).

• se tutti i coefficienti a0, a1 . . . an sono nulli, allora p(x) e detto il polinomio nullo.

Con Pn si denota l’insieme di tutti i polinomi p con grado non superiore ad n, insieme alpolinomio nullo.La principale ragione della popolarita delle funzioni polinomiali nel calcolo numerico derivadal fatto che e facile calcolarle, cioe e facile valutarle in un punto, sommarle,moltiplicarlederivarle e integrarle. La definizione mostra in particolare, che per valutare un polinomioin un punto e necessario solo moltiplicare e sommare numeri reali un numero finito divolte.

Essi giocano un ruolo centrale nella teoria dell’approssimazione e nel calcolo numericoper le loro ben note proprieta

• i polinomi di grado ≤ n formano uno spazio vettoriale di dimensione finita n + 1

• sono funzioni regolari

• sono facilmente memorizzabili, calcolabili e manipolabili su un calcolatore

• la derivata e l’antiderivata di un polinomio sono ancora polinomi i cui coefficientisono determinati algebricamente ( e anche da un calcolatore)

• il numero di zeri di un polinomio di grado n non puo superare n,

• data una qualunque funzione continua su un intervallo [a, b], esiste un polinomio chela approssima uniformemente qualunque sia l’ε prefissato: |f(x)−p(x)| < ε ∀x ∈ [a, b]

Teorema: Pn e uno spazio vettoriale di dimensione n + 1.Teorema: ( fondamentale dell’algebra) Sia p(x) un polinomio di grado n ≥ 1. L’equazione

algebrica di grado n, p(x) = 0 ha almeno una radice reale o complessa.Corollario: Ogni equazione algebrica di grado n ha esattamente n radici reali o

complesse, ciascuna contata con la sua molteplicita, cioe

p(x) = an(x− α1)m1 · (x− α2)m2 · (x− α3)m3 . . . (x− αk)mk

2011/01/12 2

dove αi, i = 1, . . . k sono le radici a due a due distinte e mi, i = 1, . . . k sono le relativemolteplicita , tali che m1 + m2 + . . .mk = n

TeoremaSiano p(x), b(x) polinomi con b(x) 6= 0; allora esistono e sono unici i polinomi q(x) ed r(x)per cui

p(x) = q(x) · b(x) + r(x)

con r(x) = 0 o r(x) con grado minore di b(x)Osservazione Dati i polinomi p(x) e b(x) con b(x) 6= 0 e sempre possibile trovare unpolinomio quoziente q(x) ed un polinomio resto r(x) per cui

p(x)b(x)

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

Valutazione numerica di un polinomioAssegnato un polinomio di grado n (nella sua forma convenzionale) p(x) =

∑ni=0 aix

i Unalgoritmo immediato per la sua valutazione in x0 e illustrato come segue: indicato con ail vettore contente i coefficienti del polinomio a = [a(n + 1), a(n), . . . , a(2), a(1)] allora

....s = 1;p = a(1);for k=2:ns = s* x_0;p = p+a(k)*s;

end

quindi il calcolo di p(x) richiede 2n noltiplicazioni e n addizioni. Se scriviamo il polinomionella forma seguente

p(x) = a0 + x(a1 + x(a2 + x(a3 + . . . + x(an−1 + xan) . . .)))

si ricava il seguente metodo dovuto a Horner:

.....p=a(n+1)for k=n:-1:1p=a(k)+x_0*p;

endr=p;

L’algoritmo di Horner richiede n moltiplicazioni e n addizioni, quindi piu rapido del prece-dente oltre che numericamente piu stabile.Lo schema di Horner e la regola di RuffiniUn caso particolarmente importante che si deduce dal teorema precedente sulla divisionedi due polinomi q(x) e b(x) e quello in cui b(x) = (x− x0) per un certo x0 reale. Allora ilrisultato dice che esistono unici i polinomi q(x) e r(x) per cui

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

dove r(x) = 0 o r(x) ha grado zero; in ogni caso r(x) e un polinomio costante diciamo r,cosı:

p(x) = q(x) · (x− x0) + r;

Se ora valutiamo p(x) in x0 si trova che p(x0) = r. Quindi vale il seguente:

2011/01/12 3

TeoremaIl resto della divisione del polinomio p(x) per (x− x0) e r = p(x0).Dal teorema del resto enunciato si deduce che per valutare un polinomio in un punto x0

si puo determinare il resto della divisione fra p(x) e (x − x0). A tal fine e nota la regoladi Ruffini:

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

x xbn xbn−1 . . . xb3 xb2 xb1

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

da cui noti i coefficienti a0, a1, . . . an ed x0 i coefficienti bi ed il resto r = p(x0) si possonoricavare cosı

bn ← an

bn−1 ← an−1 + x0bn

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

......

b1 ← a1 + x0b2

r ← a0 + x0b1

con q(x) =∑n−1

i=0 bi+1x = b1 + b2x + . . . bnxn−1; ma questo schema e esattamente quellodi Horner appena visto. Se si vogliono calcolare anche i coefficienti del polinomio q(x)l’algoritmo precedente viene cosı modificato

....p =a(n+1);b(n)=p;for k=n:-1:2

p=a(k)+p*x_0;b(k-1)=p;

endr=a(1)+x_0*b(1);

EsempioSia dato p(x) = 1 + x − 2x2 + 3x4 e si voglia calcolare p(2). Quando la valutazione sieffettua manualmente la regola di Ruffini e molto comoda:

3 0 −2 1 1

2

3

6 20 42

6 10 21 43 = P(2)

12

Valutazione di un polinomio e della sua derivataSi voglia ora determinare il valore della derivata prima di p(x) per x = x0 cioe p′(x0).Da quanto detto vale:

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

derivandop′(x) = q′(x)(x− x0) + q(x)

e per x = x0 p′(x0) = q(x0) dove q(x) =∑n−1

i=0 bi+1xi e i coefficienti bi per i = 1, 2, . . . , n

sono quelli che si ottengono applicando l’algoritmo di Ruffini Horner per valutare p(x0).Esempio

Si valuti p′(2) per p(x) = 1 + x− 2x2 + 3x4 differenziando esplicitamente p(x) si ha:

p′(x) = 1− 4x + 12x3, p′(2) = 1− 8 + 96 = 89

2011/01/12 4

3 0 −2 1 1

2

3

6 20 42

6 10 21 43 = P(2)

12

62

3 12

24

34

68

89 = P’(2)

Il calcolo di p(x0)e p′(x0) possono essere organizzati in modo da non memorizzare nessunrisultato intermedio:

....p_d=0;p=a(n+1)for k=n:-1:1p_d=p+x_0 *p_dp=a(k)+x_0*pend

.......

Analogamente si possono calcolare anche le derivate di ordine maggiore. Derivando duevolte p(x) = q(x)(x− x0) + r si ottiene

p′(x) = q′(x)(x− x0) + q(x), p”(x) = q”(x)(x− x0) + 2q′(x),

da cui p”(x0) = 2q′(x0). Quindi per calcolare p”(x0) e necessario calcolare q′(x0) utiliz-zando lo schema di Horner.

Esempio Si valuti p”(2) per p(x) = 1 + x− 2x2 + 3x4 In generale possiamo scrivere la

3 0 −2 1 1

2

3

6 20 42

6 10 21 43 = P(2)

12

62

3 12

24

34

68

89 = P’(2)

23

618

3670 = P"(2)/2

relazione fra un polinomio dividendo e il suo polinomio quoziente utilizzando un indicepj(x) = pj+1(x)(x− x0) + p(j)(x0), j = 0, . . . , n− 1, allora

• p(x) = p0(x)

• pj(x0) = j!pj(x0)

• i coefficienti di pj(x) sono i valori determinati applicando il metodo di Horner alpolinomio pj−1(x).

2011/01/12 5

Esistono varie funzioni Matlab che manipolano i polinomi:POLYVAL, POLYFIT, POLY, ROOTS, CONVSi voglia valutare p(x) = 3x7 + 2x2 + 1 nei punti equidistanti xk = −1 + k × 0.25 perk = 0, 1, . . . , 8, si puo utilizzare il seguente codice:À p = [3, 0, 0, 0, 0, 2, 0, 1]; x = [−1 : 0.25 : 1]À y = polyval(p, x)y = 0 1.7245 1.4766 1.1248 1.0000 1.1252 1.5234 2.5255 6.0000La function polyfit applica i minimi quadrati per trovare il polinomio di best-fit. Dal man-uale in linea di Matlab:

POLYFIT Fit polynomial to data.P = POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) ofdegree N that fits the data Y best in a least-squares sense. P is arow vector of length N+1 containing the polynomial coefficients indescending powers, P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).

[P,S] = POLYFIT(X,Y,N) returns the polynomial coefficients P and astructure S for use with POLYVAL to obtain error estimates forpredictions. S contains fields for the triangular factor (R) from a QRdecomposition of the Vandermonde matrix of X, the degrees of freedom(df), and the norm of the residuals (normr). If the data Y are random,an estimate of the covariance matrix of P is (Rinv*Rinv’)*normr^2/df,where Rinv is the inverse of R.

Nel caso in cui i vettori X, Y contenenti i dati da approssimare siano di dimensione Nla function polyfit calcola il polinomio di interpolazione. Come si legge dal manuale questoviene calcolato costruendo la matrice di Vandermonde che come e noto e mal condizionata,per cui se il set di dati genera una matrice V malcondizionata e molto probabile ottenereil seguente messagio:

Warning: Polynomial is badly conditioned. Add points with distinct Xvalues, reduce the degree of the polynomial, or try centeringand scaling as described in HELP POLYFIT.

Matlab fornisce i coefficienti del polinomio calcolato, ma informa che poiche’ la matricedei coefficienti e mal condizionata i dati possono non essere attendibili.

Con il comando p = conv(p1, p2) si calcolano i coefficienti del polinomio ottenutocome prodotto dei polinomi i cui coefficienti sono precisati in p1 e p2. Invece il comando[q, r] = deconv(p1, p2) calcola i coefficienti del quoziente e del resto della divisione fra p1e p2, cioe q ed r tali chep1 = conv(p2, q) + r

EsempioConsideriamo i polinomi p1(x) = x4 − 1 e p2(x) = x3 − 1, calcoliamo il prodotto e ladivisioneÀ p1 = [1, 0, 0, 0,−1];À p2 = [1, 0, 0,−1];À p = conv(p1, p2)p =1, 0, 0,−1,−1, 0, 0, 1À [q, r] = deconv(p1, p2)

2011/01/12 6

q =1 0r =0 0 0 1 − 1Troviamo pertanto i polinomi p(x) = p1(x)p2(x) = x7−x4−x3+1, q(x) = x e r(x) = x−1tali che p1(x) = q(x) ∗ p2(x) + r(x).

POLY(A), when A is an N by N matrix, is a row vector withN+1 elements which are the coefficients of thecharacteristic polynomial, DET(lambda*EYE(SIZE(A)) - A) .

Esempio

A =

1/3 2/3 3 01 0 −1 2 00 0 −5/3 −2/3 00 0 1 0 1/20 0 0 0 2/3

À p = poly(A)p = 1 2/3 − 13/9 − 962/999 4/9 296/999]quindi il polinomio p(x) ha la seguente espressione:p(x) = x5 + 2/3x4 − 13/9x3 − 962/999x2 + 4/9x + 296/999

POLY(V), when V is a vector, is a vector whose elements arethe coefficients of the polynomial whose roots are theelements of V . For vectors, ROOTS and POLY are inversefunctions of each other, up to ordering, scaling, androundoff error.

À c = poly([1, 1, 1, 1, 2])c = 1 − 6 14 − 16 9 − 2quindi il polinomio ha la seguente espressione:p(x) = x5 − 6x4 + 14x3 − 16x2 + 9x− 2

ROOTS Find polynomial roots.ROOTS(C) computes the roots of the polynomial whose coefficientsare the elements of the vector C. If C has N+1 components,the polynomial is C(1)*X^N + ... + C(N)*X + C(N+1).

Note: Leading zeros in C are discarded first. Then, leadingrelative zeros are removed as well. That is, if division by theleading coefficient results in overflow, all coefficients up to thefirst coefficient where overflow ocurred are also discarded. Thisprocess is repeated until the leading coefficient is’nt a relative 0

À r = roots(c)

r= 2.00001.0002 + 0.0002i1.0002 - 0.0002i0.9998 + 0.0002i0.9998 - 0.0002i

>> S=compan(c)S=

2011/01/12 7

6 -14 16 -9 21 0 0 0 00 1 0 0 00 0 1 0 00 0 0 1 0

>> eig(S)ans = 2.0000

1.0002 + 0.0002i1.0002 - 0.0002i0.9998 + 0.0002i0.9998 - 0.0002i

2011/01/12 8

Interpolazione

Di una funzione f(x) siano noti i valori fi = f(xi) per i = 0, 1, . . . , n, in corrispondenzadi n + 1 punti distinti xi. Si vuole determinare una funzione φ tale che

φ(xi) = fi i = 0, 1, . . . , n

Una tale funzione φ si dice che interpola la f nei punti xi i = 0, 1, . . . , n.Nel caso generale viene individuato uno spazio di dimensione finita P di funzioni realie una loro base di rappresentazione ϕj(x), j = 0, 1, . . . , n per cui ogni elemento di Prisulta definito da una combinazione lineare a coefficienti reali delle ϕj :

n∑

j=0

cjϕj(x)

allora il problema di interpolare i dati assegnati con una funzione φ ∈ P consiste neldeterminare la funzione φ∗ o meglio i coefficienti c∗0, c∗1, . . . c∗n che la definiscono tale che

φ∗(xi) = fi i = 0, 1, . . . , n

La scelta dello spazio solitamente e dettata dal fatto che esista sempre una soluzioneunica al problema, successivamente dalla regolarita desiderata per l’interpolante e dallaparticolarita dei dati.

Abbiamo visto che risulta molto importante scegliere la classe delle funzioni ϕj(x) siaper quanto riguarda il condizionamento sia per quanto riguarda il modello dell’approssimazione.La scelta deve essere motivata in modo da tener conto delle specifiche proprieta della f(x).Inoltre devono essere facilmente calcolabili e dotate di buone proprieta di regolarita. Gen-eralmente si utilizzano:

• la classe delle funzioni razionali

• la classe delle funzioni trigonometriche.

• la classe delle funzioni esponenziali

Occorre tener presente che non sempre e consigliabile determinare una funzione approssi-mante la f richiedendo che abbia gli stessi valori in certi punti. Se infatti della f sononoti solo approssimazioni dei valori della f , per esempio i dati provengono da misurazioni,e meglio usare un metodo di approssimazione che possa attenuare gli errori dovuti allemisurazioni.In queste lezioni tratteremo solamente il caso dell’interpolazione polinomiale che garan-tisce sempre l’esistenza e l’unicita della soluzione, una buona regolarita ed una buonaricostruzione in casi particolari.

TeoremaDati (n + 1) punti di interpolazione (xi, fi), i = 0, 1, . . . , n con xj 6= xk k 6= j esiste ed eunico il polinomio p ∈ Pn che verifica le condizioni

p(xi) = f(xi)

DimostrazioneSi consideri il vettore a = (a0, a1, . . . , an)T dei coefficienti del polinomio p(x) = a0 +a1x+a2x

2 + . . . anxn, il vettore f = (f0, f1, . . . , fn)T e la matrice V = {νi,j} detta matrice diVandermonde, cosı definita:

νi,j = xj−1i−1 i = 1, 2, . . . , n + 1 j = 1, 2, . . . , n + 1

2011/01/12 9

Imponendo che p(x) verifichi le n + 1 condizioni di interpolazione, si ottiene il sistemalineare di n + 1 equazioni in n + 1 incognite

V a = f.

La matrice di vandermonde risulta non singolare in quanto

detV = Πni,j=0;j>i(xj − xi)

e i punti xi sono per ipotesi distinti, ne segue che il sistema lineare ha una ed una solasoluzione e quindi il polinomio p(x) esiste ed e unico.

Osservazione Solo raramenmte del polinomio di interpolazione sono richiesti i coef-ficienti, in generale si vuole il valore di p(x) in uno o piu punti, e per calcolarlo non econveniente risolvere il sistema lineare, perche cio richiederebbe un numero di operazionidell’ordine di n3/3 operazioni moltiplicative, inoltre la matrice di Vandermonde puo esseremal condizionata e quindi il calcolo di p(x) sarebbe numericamente instabile.Il polinomio di interpolazione pur essendo unico puo essere rappresentato in diverse formepiu convenienti sia dal punto di vista del costo computazionale che della stabilita.

Base di Lagrange

Dati (n + 1) punti di interpolazione (xi, yi), i = 0, 1, . . . , n con a = x0 < x1 < . . . <xn = b, il polinomio di interpolazione si scrive nella forma di Lagrange come segue:

p(x) =n∑

i=0

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

Le funzioni di Lagrange sono definite per i = 0, 1, . . . , n:

Li,n(x) =x− x0

xi − x0× . . .× x− xi−1

xi − xi−1× x− xi+1

xi − xi+1× . . .× x− xn

xi − xn

= Πnj 6=i,j=0

x− xj

xi − xj

Le funzioni di Lagrange Li,n(x) formano una base per lo spazio Pn dei polinomi digrado n

{L0,n(x), L1,n(x), . . . Ln,n(x)}Per come sono definite valgono le seguenti proprieta

Li,n(xj) = δij =

{1 se j = i0 se j 6= i

Li,n(x) e un polinomio di grado al piu n, con n zeri distinti: x0, x1, . . . xi−1, xi+1, . . . xn.Quindi determinare il polinomio p(x) che soddisfi le condizioni di interpolazione p(xi) = yi,∀i, nella forma di Lagrange

p(x) = a0L0,n(x) + a1L1,n(x) + . . . + anLn,n(x)

imponendo le condizioni di interpolazione si ottiene il sistema lineare

L0,n(x0) L1,n(x0) . . . Ln,n(x0)L0,n(x1) L1,n(x1) . . . Ln,n(x1)...

......

...L0,n(xn) L1,n(xn) . . . Ln,n(xn)

·

a0

a1...an

=

y0

y1...yn

2011/01/12 10

Si osservi che la matrice L risulta uguale alla matrice identita di ordine n + 1, In+1.Quindi nell’interpolazione secondo Lagrange non cerchiamo i coefficienti ai perche sonoidenticamente uguali ai valori noti yi. In questo caso cerchiamo proprio le funzioni base diLagrange costruite sulle ascisse xi. L’idea e di risolvere gli n+1 problemi di interpolazionedove le ascisse x0, x1, . . . , xn sono quelle del problema iniziale in corrispondenza delle qualisi hanno n + 1 insiemi di valori yi cosı fatti:{1, 0, 0, . . . , 0} ⇒ unica soluzione L0,n(xj){0, 1, 0, . . . , 0} ⇒ unica soluzione L1,n(xj){0, 0, 1, . . . , 0} ⇒ unica soluzione L2,n(xj)...{0, 0, 0, . . . , 1} ⇒ unica soluzione Ln,n(xj)La soluzione unica al problema iniziale si ottiene combinando assieme le n + 1 soluzioniparziali:

p(x) =n∑

i=0

yiLi,n(x)

e il polinomio unico di interpolazione e

Li,n(xj) = δi,j ⇒ p(xi) =n∑

i=0

yiLi,n(xi) = yi ∀i.

EsempioDati i punti x0 = 0, x1 = 1, x2 = 3 Le funzioni base di Lagrange ( anche dette funzionicardinali) costruite su queste ascisse sono

L0,2 =x− 10− 1

· x− 30− 3

=(x− 1)(x− 3)

3=

x2 − 4x + 33

L1,2 =x− 01− 0

· x− 31− 3

=−x(x− 3)

2

L2,2 =x− 03− 0

· x− 13− 1

=x(x− 1)

6Il polinomio che interpola i tre punti (0, y0), (1, y1), (3, y2) e allora

0 0.5 1 1.5 2 2.5 3−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

p(x) = y0x2 − 4x + 3

3− y1

−x(x− 3)2

+ y2x(x− 1)

6

Le funzioni base di Lagrange costruite su n+1 ascisse distinte hanno grado≤ n nell’esempioprecedente le tre funzioni base di Lagrange hanno tutte grado 2. Se i punti (xi, yi) fossero

2011/01/12 11

tutti su una retta, la combinazione lineare si ridurrebbe in modo che il polinomio p risul-tante abbia grado 1.EsempioDati i punti (0, 0), (1, 1), (3, 3), le funzioni base di Lagrange sono:

L0,2 =x2 − 4x + 3

3, L1,2 =

−x(x− 3)2

i e L2,2 =x(x− 1)

6

p(x) = 0x2 − 4x + 3

3− x(x− 3)

2+ 3

x(x− 1)6

= x.

Se si deve valutare il polinomio interpolatore espresso nella forma di Lagrange p(x) =∑i yiLi,n(x) in un punto xi 6= xi, ∀i, si e tentati dal volerlo trasformare nella forma mono-

miale p(x) =∑

i aixi per poi applicare l’algoritmo di Ruffini-Horner. Procedere in questo

modo, vuole dire dover calcolare gli ai che richiedono oltre 6n2 operazioni (×, +)Un altro algoritmo potrebbe essere:Dati x, n, {x0, x1, . . . , xn}, {y0, y1, . . . , yn} p ← 0for j=0,1,. . .,n....... q ← 1....... for i=0,1,. . .,n...............if j 6= i............... q ← q × (x−xi)

(xj−xi)

....... end

....... p ← p + yj × qendma anche questo ha un costo di circa 4n2 operazioni (×, +).

Per cui occorre cambiare strategia.Consideriamo

ω(x) = (x− x0)(x− x1) . . . (x− xn−1)(x− xn) = Πnk=0(x− xk)

Quindi possiamo scrivere

p(x) =n∑

i=0

yiLi,n(x) = ω(x)n∑

i=0

zi

x− xi

conzi =

yi

Πnj 6=i,j=0(xi − xj)

Infattin∑

i=0

yiLi,n(x) =n∑

i=0

yiΠnj 6=i,j=0

x− xj

xi − xj=

=n∑

i=0

yi

Πnj 6=i,j=0(x− xj)

Πnj 6=i,j=0(xi − xj)

=n∑

i=0

yiω(x)/(x− xi)

Πnj 6=i,j=0(xi − xj)

=

=n∑

i=0

yi

Πnj 6=i,j=0(xi − xj)

ω(x)x− xi

=n∑

i=0

ziω(x)

x− xi

2011/01/12 12

= ω(x)n∑

i=0

zi

x− xi

ove abbiamo postozi =

yi

Πnj 6=i,j=0(xi − xj)

Quindi per calcolare

p(x) =n∑

i=0

yiLi,n(x) = ω(x)n∑

i=0

zi

x− xi

Le operazioni richieste sono, a meno di termini di ordine inferiore,

• n sottrazioni per calcolare (x− xi), per i = 0, 1, . . . , n

• n2/2 sottrazioni per xj − xi, per i = 0, 1, . . . , n, e j = 0, 1, . . . , n,

• n2 prodotti per calcolare zi, per i = 0, 1, . . . , n;

• n somme e 2n prodotti per calcolare p(x).

In totale questo calcolo richiede n2/2 addizioni ed n2 prodotti. Per calcolare il valoredello stesso polinomio p in un altro punto x conviene sempre usare questa forma perchegli scalari zi dipendono solo da xi e possono essere calcolati una volta per tutte. Quindi ilcalcolo in un altro punto x richiede solo il passo finale cioe n somme e 2n prodotti. Infinelasciando inalterati gli xi e variando solo le ordinate yi si puo ancora calcolare p(x) con nsomme e 2n prodotti aggiuntivi.

2011/01/12 13

RIEPILOGOI costi computazionali dell’interpolazione con la forma di Lagrange sono i seguenti:Calcolo in un punto: n2/2 somme ed n2 prodottiCalcolo in m punti:n2/2 + 2nm somme ed n2 + 2mn prodotti.In MATLABy = lagrval(y, x, x)

function yy=lagrval(x,y,xx)% valutazione di un polinomio nella forma di Lagrange% y --> vettore dei coefficienti% x --> vettore dei punti in base ai quali restano determinati% i polinomi elementari di Lagrange% xx --> vettore di ascisse in corrispondenza delle quali% si vuole valutare il polinomio% yy <-- vettore dei valori del polinomion=length(x);den(1)=prod(x(1)-x(2:n));for i=2:n-1den(i)=prod(x(i)-x(1:i-1))*prod(x(i)-x(i+1:n));

endden(n)=prod(x(n)-x(1:n-1));m=length(xx);for i=1:momega=prod(xx(i)-x);if omega==0

[nz k]=max(xx(i)==x);yy(i)=y(k);

elseyy(i)=omega*sum(y./(den.*(xx(i)-x)));

endend

P. di Lagrange P.di Newton1 punto n2/2(+) + n2(×) n2(+) + n2/2(×)

in m punti n2/2 + 2mn)(+)+ (n2 + 2mn)(+)++(n2 + 2mn)(×) +n2/2 + 2mn)(×)

2011/01/12 14

Analisi dell’errore di interpolazione

Dati n+1 osservazioni fi di una funzione f(x) in corrispondenza di n+1 ascisse distinteabbiamo visto alcuni metodi per costruire il polinomio di grado minimo che interpola ipunti (xi, fi). In questo caso ha senso chiedersi che errore si commette in un puntox nel considerare il valore di p(x) al posto di f(x), cioe quanto e grande l’errore diinterpolazione

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

che si commette in un puntox ∈ [a, b], x 6= xi, ∀i.

Per dare una risposta e necessario porre alcune ipotesi di regolarita sulla f .

TeoremaSia f(x) ∈ Cn+1

[a,b] ossia f continua con tutte le sue derivate fino a quella di ordine n + 1.Sia [a, b] un intervallo di R limitato e chiuso contenente le ascisse distinte xi. Sia x ∈ [a, b]un punto qualsiasi, distinto dagli xi. Esiste allora un punto ξ ( dipendente da x) internoad [a, b] per cui:

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

(n + 1)!fn+1(ξ)

conω(x) = (x− x0)(x− x1) . . . (x− xn) = Πn

k=0(x− xk).

CorollarioSi indichi con

Mn+1 = maxx∈[a,b]

|fn+1(x)|

allora un limite superiore ad R(x) = f(x)− p(x) e dato da

|R(x)| = |ω(x)|(n + 1)!

Mn+1.

Esempiof(x) = ex, f ∈ Cn+1

[a,b] , si ha:

Mn+1 = maxx∈[a,b]

|ex| = eb

Inoltre per ogni scelta delle ascisse xi ∈ [a, b] vale|ω(x)| = |(x− x0)(x− x1) . . . (x− xn−1)(x− xn)|≤ (b− a)(b− a) . . . (b− a) = (b− a)n+1 Pertanto

maxx∈[a,b]

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

(n + 1)!eb

ma in questo caso

limn→∞ max

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

n→∞(b− a)n+1

(n + 1)!eb = 0

quindi in altre parolelim

n→∞ |p(x)− f(x)| = 0

Il polinomio p di interpolazione converge uniformemente ad f su [a, b] quando il numeron di punti xi di interpolazione tende a ∞.

2011/01/12 15

Estrapolazione L’espressione dell’errore R(x) evidenzia come nel caso di estrapo-lazione

x /∈ [a, b] := [mini

xi, maxi

xi]

l’errore R(x) puo aumentare notevolmente, assieme a |ω(x)|.

Errori nell’interpolazione, limitazioni dell’errore, convergenza e punti diChebyshevNella trattazione del problema di interpolazione sino ad ora abbiamo trattato le ordinateyi dei punti di interpolazione (xi, yi) come numeri arbitrari. Assumiamo che soddisfinola relazione f(xi) = yi dove f abbia tutte le derivate necessarie. Quindi si dice che ilpolinomio p(x) calcolato interpola la f in x0, x1, . . . , xn. Tale polinomio viene spessocostruito per trovare una approssimazione della f piu semplice da trattare e naturaleconsiderare la seguente espressione dell’errore:

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

in seguito assumiamo che la variabile x non assuma valori uguali a x0, x1, . . . , xn ( perchein questi punti l’errore e zero).

Aumentando il numero di punti di interpolazione aumenta in corrispondenza il gradodel polinomio, al crescere quindi di n la successione

{pn(x)}

puo non convergere ad f(x) puo capitare anche nel caso in cui f ∈ C∞[a,b]

Consideriamo la seguente funzione di Runge

f(t) = 1/(1 + t2)

se interpolata con punti equispaziati nell’intervallo [−5, 5] l’interpolante diverge in cor-rispondenza degli estremi dell’intervallo. Si puo dimostrare che al crescere di n la succes-sione dei polinomi di interpolazione {p(n} costruiti sulle ascisse equidistanti:

xi = i10n− 1, i = 0, 1, . . . , n

non converge alla funzione di Runge, perche gli errori diventano arbitrariamente grandispecialmente nei punti vicino agli estremi dell’intervallo. Quindi imporre la regolarita sullaf non e sufficiente a garantire la convergenza dell’interpolazione.Se si considerano come ascisse di interpolazione gli n + 1 zeri del polinomio di Chebyshevdi grado n rispettivamente:xi = cos( (2i+1)π

2(n+1) ) ∈ [−1, 1] i = 0, 1, . . . , noppurexi = a+b

2 − b−a2 cos( (2i+1)π

2(n+1) ) ∈ [a, b] i = 0, 1, . . . , nCon questa distribuzione delle ascisse xi si dimostra che, al crescere di n →∞, il polinomiopn(x) di interpolazione converge a f(x) quando x ∈ [a, b].

2011/01/12 16

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

0

0.5

1interpolazione funzione di Runge punti equidistanti

grado del polinomio 6−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.5

0

0.5

1

1.5

2interpolazione funzione di Runge punti equidistanti

grado del polinomio 10

Questo comportamento e dovuto al polinomio ω(x) = (x− x0)(x− x1) . . . (x− xn) checompare nell’espressione dell’errore che assume valori di picco agli estremi dell’intervallo.L’errore di interpolazione puo essere limitato se si scelgono i punti di interpolazione inmodo che

| maxx∈[a,b]

ω(x)|

sia minimo. I punti di Chebyshev sono un tentativo di aggiustare i punti di interpolazione

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

0

0.5

1interpolazione funzione di runge

grado del polinomio 6−5 −4 −3 −2 −1 0 1 2 3 4 5

−0.5

0

0.5

1interpolazione funzione di runge

grado del polinomio 10

poli. di grado 6 interpolante nei punti di Chebyshev evidenziati poli. di grado 10 interpolante nei punti di Chebyshev evidenziati

Si puo dimostrare cheminω(x)

maxx∈[a,b]

|ω(x)| = (b− a)−n

e il minimo e raggiunto quando

xi =a + b

2− b− a

2cos(

(2i + 1)π2(n + 1)

) ∈ [a, b] i = 0, 1, . . . , n

che sono i punti di Chebyshev.

2011/01/12 17

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

−3000

−2000

−1000

0

1000

2000

3000

4000grafico della funzione omega

grado del polinomio 6

equi.chebychev

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

−4

−3

−2

−1

0

1

2

3

4

5x 10

5 grafico della funzione omega

grado del polinomio 10

equi.chebychev

funzione ω per polinomi di grado 6 funzione ω per polinomi di grado 10

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

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1interpolazione funzione di runge

grado del polinomio 20−5 −4 −3 −2 −1 0 1 2 3 4 5

−1.5

−1

−0.5

0

0.5

1

1.5x 10

11 grafico della funzione omega

grado del polinomio 20

equi.chebychev

Poli. di grado 20 interpolante nei punti di Chebyshev evidenziati funzione ω per polinomi di grado 20

Sfortunatamente ci sono funzioni per cui l’interpolante nei punti di Chebyshev nonconverge e occorre cambiare strategia.

2011/01/12 18

Interpolazione polinomiale a tratti

Si intende l’interpolazione di un set di dati su un intervallo con piu polinomi ciascunodei quali definito in un sottointervallo dell’intervallo dato. In particolare siano assegnatem+1 osservazioni in corrispondenza di m+1 punti distinti e ordinati, cioe siano assegnati(xi, yi), i = 0, . . . , m con xi < xi+1. Allora una interpolante polinomiale a tratti consistein m polinomi pi(x), i = 0, . . . , m− 1 di grado n, definiti sugli intervalli [xi, xi+1] con leseguenti 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 una funzione continuafino alla derivata di ordine k.

Polinomi a tratti

Abbiamo osservato che i polinomi hanno un comportamento inaccettabile quando ven-gono utilizzati per approssimare dati o funzioni. Quando superano il quinto o sesto gradooscillano fortemente su intervalli di definizione grandi, mentre hanno un buon comporta-mento se di grado basso e su intervalli piccoli. Questo fenomeno suggerisce di utilizzarepolinomi di grado relativamente basso e suddividere l’intervallo in piccole parti. Si puodare la seguente definizione:

Definizione: Polinomio a TrattiSia [a, b] un intervallo limitato e chiuso e sia ∆ la partizione dell’intervallo [a, b] data da:

∆ = {xi}i=0,...,m

su un insieme di punti , detti nodi, tali che:

a = x0 < x1 < . . . < xm = b

La partizione di [a, b] indotta dall’insieme ∆ risulta

• I0 = [x0, x1]

• Ii = [xi, xi+1] i = 1, . . . , m− 1

• Im = [xm−1, xm]

Indichiamo con Pn lo spazio dei polinomi a coefficienti reali di ordine al piu n, definiamolo spazio dei polinomi a tratti:

T Pn = {ϕ | ∃p0, p1, . . . pm ∈ Pn; tale cheϕ(x) = pi(x),

∀x ∈ Ii, i = 0, 1, . . . , m}Esempio

Consideriamo il polinomio a tratti di grado 1 che interpola la funzione di Runge f(x) =1

1+x2 nei 7 nodi equidistanti assegnati in figuraSi osserva che utilizzando i polinomi a tratti si e guadagnato in flessibilita, ma si e perso

in regolarita: l’interpolante non ha piu’ la derivata prima continua. In molte applicazionisi richiede che la funzione approssimante sia almeno di classe C1.

2011/01/12 19

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

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1interpolazione a tratti con polinomi lineari

Interpolazione locale

Per interpolazione locale si intende un metodo che individua il polinomio interpolantea tratti determinando ogni singolo polinomio pi(x) in modo indipendente dagli altri.Consideriamo il caso particolare di polinomi a tratti cubici, con regolarita C1; i singolipolinomi cubici pi(x) possono essere espressi nella base di Bernstein e determinanti inmodo da interpolare le ordinate yi e yi+1 in corrispondenza di xi e xi+1 sull’intervallo[xi, xi+1] e i valori di derivata y′i e y′i+1 sempre in corrispondenza di xi e xi+1. In questocaso si parla di interpolanti cubiche di Hermite. I metodi di interpolazione locale fannouso delle derivate y′i in corrispondenza dei punti assegnati xi. Se queste non sono date,devono essere calcolate e diventano parte integrante del problema di interpolazione.

Interpolazione globale

Per interpolazione globale si intende un metodo che determina il polinomio a trattiinterpolante ricavando ciascun polinomio pi(x) in modo dipendente dagli altri. Nel casoparticolare di polinomi a tratti cubici con massima regolarita C2 si parla di Spline cu-biche di Interpolazione.

Lo spazio delle funzioni Spline polinomiali di grado n puo essere sinteticamente scrittocome :

Sn(∆) = T Pn(∆) ∩ C(n−1)[a,b]

Questa condizione implica che ogni elemento di Sm(∆) e sufficientemente flessibile, perchesi hanno funzioni continue e derivabili con continuita sino all’ordine n−1 che e la massimaregolarita che si puo richiedere ai polinomi a tratti di grado m-1.Consideriamo le funzioni Spline cubiche.Riassumendo:siano x0 < x1 < x2 < . . . < xm punti dati detti nodi , una Spline Cubica e una funzioneS(x), definita sull’intervallo [x0, xm] con le seguenti proprieta:

• e una funzione con derivate prime e seconde continue

• ristretta ad ogni intervallo Ii e un polinomio di grado 3

• S(xj) = f(xj), per j = 0, 1, 2, . . . , m

2011/01/12 20

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

S0

S1

Sm−1S

m−2

x0

x1

x2

xm−1 x

m

Se poniamo S(x) = si(x) per xi ≤ x ≤ xi+1, gli si(x) saranno m polinomi cubici deltipo si(x) = ai + bix + cix

2 + dix3 per x ∈ [xi, xi+1], i = 0, 1, . . . , k per ciascuna cubica

occorrono 4 condizioni quindi in totale 4m condizioni per determinare tutti i coefficienti.Ogni polinomio interpola la f negli estremi del proprio intervallo. Questo da 2m condizioni.Le condizioni di continuita delle drivate nei nodi interni rispettivamente:

• s′i−1(xi) = s′i(xi) per i = 1, 2, . . . ,m− 1

• s”i−1(xi) = s”i(xi)

sono altre 2(m− 1) condizioni.In tutto si hanno 2m + 2(m− 1) = 4m− 2 condizioni per detereminare i 4m coefficienti.Le restanti due condizioni vengono scelte fra:

1. s′0(x0) = y′0 e s′m−1(xm) = y′m, con y′0 e y′m assegnati dall’utente Spline con derivateagli estremi,

2. s”(x0) = s”(xm) = 0 Spline di Interpolazione Naturali

3. s0(x0) = sm(xm), s′(x0) = s′(xm) e s”(x0) = s”(xm) Spline di Interpolazione peri-odiche (si usano quando y0 = ym, cioe’ quando si vogliono interpolare curve chiuseo periodiche).

In alternativa alla seconda condizione, si puo richiedere che la derivata terza di S(x) neinodi x0 e xm sia continua. Questa condizione a cui viene attribuito il nome not-a-knot-condition e quella implementata nel comando Matlab spline ( si veda anche il toolboxSplines).Per determinare i coefficienti ai, bi, ci, di si potrebbe al solito usare tutte le condizionielencate per ricavare un sistema avente come incognite i coefficienti e risolverlo, ma non eil modo piu efficiente. Un metodo efficiente si ottiene prendendo come incognite i momenticioe le derivate seconde nei nodi xi

Mi = S”3(xi)

essendo ogni S”(x) una spline del primo ordine cioe una funzione lineare a tratti avremo

S”(x) =(x− xi−1)Mi + (xi − x)Mi−1

xi − xi−1x ∈ [xi−1, xi]

Posto hi = (xi − xi−1) integrando due volte nell’intervallo [xi−1, xi]si ha

S′(x) =(x− xi−1)2Mi − (xi − x)2Mi−1

2hi+ Ci x ∈ [xi−1, xi]

2011/01/12 21

e

S(x) =(x− xi−1)3Mi + (xi − x)3Mi−1

6hi+ Ci(x− xi−1) + Di x ∈ [xi−1, xi]

imponendo le condizioni di interpolazione si(xi) = si−1(xi) si ottengono i parametri Ci eDi in funzione dei momenti Mi piu precisamente

(xi−xi−1)3Mi−1

6hi+ Di = yi− 1 i = 1, . . . , m

(xi−xi−1)3Mi

6hi+ Ci(xi − xi−1) + Di = yi

da cui si ricava

Di = yi−1 − h2i

6Mi−1, Ci =

(yi − yi−1)hi

− hi(Mi −Mi−1)

6

infine sostituendo nella S′(x) si ottiene per x ∈ [xi−1, xi]

S′(x) =(x− xi−1)3Mi + (xi − x)3Mi−1

6hi+

yi − yi−1

2hi− hi

Mi −Mi−1

6

Imponendo la continuita di S′(x) nei nodi (rispettivamente s′i(xi) = s′i−1(xi)) per i =1, . . . ,m− 1 si ha

hiMi−1 + 2(hi + hi+1)Mi + hi+1Mi+1 = 6(

yi+1 − yi

hi+1− yi − yi−1

hi

).

Tenendo conto che per le spline naturali M0 = Mm = 0 si ottiene un sistema di equazionilineari nelle m− 1 incognite M1, M2, . . .Mi−1 la cui matrice dei coefficienti e

2(h1 + h2) h2 . . . 0h2 2(h2 + h3) . . . 0

0. . . . . .

......

. . . hm−1

0 . . . hm−1 2(hm−1 + hm)

La risoluzione di questo sistema comporta un costo computazionale dell’ordine di 3n molti-plicazioni. Si osservi che nel caso in 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

D’ora in poi per calcolare una spline utilizzeremo il comando MATLABspline(x,y,z) i parametri di ingresso di questa function sono i vettori x e y dei nodi diinterpolazione ed il vettore z delle ascisse nelle quali si vuole che venga valutata la spline.

Le curve spline sono usate in numerose applicazioni industriali e sono alla base deimoderni sistemi CAD, il vantaggio di usare una spline e di fornire movimenti regolarievitando bruschi cambi di velocita e accelerazione se regolare sino alla derivata seconda.

Riprendiamo l’esercizio riguardante l’interpolazione della funzione di Runge:

f(x) =1

x2 + 1

.

2011/01/12 22

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

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Interpolazione mediante spline cubica naturale della funzione Runge con 19 punti di interpolazione

Funzione RungeSpline cubica interpolantePunti di interpolazione

Figure 1: Interpolazione della funzione di Runge con una spline cubica naturale su una partizione contenente 19 punti

nell’intervallo [−5; 5].

Interpolazione con il polinomio di grado 6 sui 7 punti equidistanti della funzione diRunge f(x) = 1

1+x2 definita sull’intervallo [−5, 5]

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

0

0.5

1interpolazione funzione di Runge punti equidistanti

grado del polinomio 6 −5 −4 −3 −2 −1 0 1 2 3 4 50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Interpolazione mediante spline cubica naturale della funzione Runge con 7 punti di interpolazione

Funzione RungeSpline cubica interpolantePunti di interpolazione

. Interpolazione della funzione di Runge con una spline cubica naturale su una par-tizione nodale di 7 punti equidistanti nell’intervallo [−5; 5].

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

0

0.5

1

1.5

2interpolazione funzione di Runge punti equidistanti

grado del polinomio 10 −5 −4 −3 −2 −1 0 1 2 3 4 50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Interpolazione mediante spline cubica naturale della funzione Runge con 11 punti di interpolazione

Funzione RungeSpline cubica interpolantePunti di interpolazione

[Interpolazione della funzione di Runge con una spline su 11 punti]Interpolazione della funzione

di Runge con una spline cubica naturale su una partizione nodale di 11 punti nell’intervallo [−5; 5].

2011/01/12 23

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

In tabella e riportato l’errore massimo per l’interpolazione polinomiale cubica a trattiC2 della funzione di Runge su punti equidistanti.

Esempio

Nella tabella seguente sono riportate le variazioni di temperatura media annua delglobo al variare della concentrazione K di acido carbonico e della latitudine:

latitudine K=0.67 K=1.5 K=2.0 K=3.065 -3.1 3.52 6.05 9.355 -3.22 3.62 6.02 9.345 -3.3 3.65 5.92 9.1735 -3.32 3.52 5.7 8.8225 -3.17 3.47 5.3 8.115 -3.07 3.25 5.02 7.525 -3.02 3.15 4.95 7.3-5 -3.02 3.15 4.97 7.35-15 -3.12 3.2 5.07 7.62-25 -3.2 3.27 5.35 8.22-35 -3.35 3.52 5.62 8.8-45 -3.37 3.7 5.95 9.25-55 -3.25 3.7 6.1 9.5

Le tecniche di approssimazione previste consentono di ricostruire a partire dai datidisponibili i valori della temperatura anche per latitudini o concentrazioni non previstedalla tabella stessa. Consideriamo la colonna K = 0.67 e calcoliamo la spline cubicainterpolatoria S3 . Per ottenere le valutazioni di S3 in corrispondenza dei valori zi =−55 + i, i = 0, . . . , 120 occorre dare le seguenti istruzioni:

À x = [−55 : 10 : 65];

À y = [−3.25,−3.37,−3.35,−3.2,−3.12,−3.02,−3.02, ...

−3.07,−3.17,−3.32,−3.3,−3.22,−3.1]

À z = [−55 : 1 : 65];

s = spline(x, y, z)

Si noti l’uso dei tre punti di continuazione tramite i quali si informa Matlab che la riga con-tenete l’istruzione continua. La linea tratteggiata e il grafico del polinomio interpolatoredi grado 12.

Errore di interpolazione

La teoria dell’errore nell’approssimazione spline e piu complessa di quella relativaall’ordinaria interpolazione a tratti, qui illustriamo i risultati principali:TeoremaFra tutte le funzioni g(x) ∈ C2

[a,b] tali che g(xi) = f(xi) per i = 0, . . . , n (cioe soddisfanole condizioni di interpolazione), la spline naturale s(x) e quella che minimizza l’integrale

∫ b

a[g”(x)]2 dx

L’integrale puo essere assunto come una misura della curvatura globale della funzione g(x),quindi risulta che la spline cubica naturale e quella che minimizza la curvatura globale. Siricorda che la curvatura descritta dall’equazione

y = f(x)

2011/01/12 24

−60 −40 −20 0 20 40 60 80−3.5

−3.4

−3.3

−3.2

−3.1

−3

−2.9

−2.8

−2.7pol.12 grad0nodispline cubica

e data da:|f”(x)| · {1 + f ′(x)2}−3/2

quindi se f ′(x) e sufficientemente piccolo, f”(x) risulta una approssimazione della cur-vatura e l’operazione di integrazione risulta una misura di quanto oscilla la funzione.Quando non si conosce nulla circa la f(x), la scelta della spline not a knot condition im-plementata sotto Matlab e la migliore. Il nome indica che il primo e l’ultimo nodo internox1, xm−1 non si comportano come nodi normali infatti se richiediamo la continuita dellas3(x) in quei punti troviamo che i primi due tratti cubici sono identici analogamente pergli ultimi due, le due condizioni di continuita equivalgono ad eliminare i nodi x2 e xm−1

dalla suddivisione mantenedo comunque le condizioni di interpolazione.Per tale motivotali consizioni sono note come not a knot condition e fornisce una esemplificazione delfatto che nel caso di funzioni polinomiali a tratti i punti della suddivisione non coincidononecessariamente con i nodi dell’interpolazione.

Dal teorema sopra esposto segue anche che se f(x) ∈ C2[a,b] allora

∫ b

a[S”(x)]2 dx ≤

∫ b

a[f”(x)]2 dx

a differenza di quanto accade per i polinomi che interpolano la f(x) su tutto l’intervallo[a, b] quando si infittiscono i nodi xi, i = 0, . . . , n le funzioni spline convergono alla f(x). Nell’ipotesi in cui la f(x) abbia derivate continue sino a quella del quart’ordine vale lasequente disuguaglianza

maxx∈I

|f r(x)− sr3(x)| ≤ CrH

4−r maxx∈I

|f4(x)|, r = 0, 1, 2, 3

ove H e la masima ampiezza degli intervalli Ii e le costanti Cr sono indipendenti da H.E quindi evidente che non solo la funzione f , ma anche le sue derivate prima, seconda eterza vengono bene approssimata dalla S(x).

2011/01/12 25

Integrazione numerica

L’integrazione numerica o quadratura numerica consiste nel calcolare un’approssimazioneper ∫ b

af(x)dx

Il problema sorge quando l’integrazione non puo essere eseguita esattamente, cioe non sipuo trovare una funzione primitiva come per esempio per

∫ 1

0e−x2

dx, o

∫ π2

0

√(1 + cos2(x))dx.

Anche nel caso in cui la si conosca, potrebbe essere complicato valutarla. Ad esempio

1π

∫ π

0cos(4π)cos(3sin(x))dx = (

32)4

∞∑

k=0

(−9/4)4

k!(k + 4)!,

in questo ultimo caso il problema del calcolo dell’integrale si e trasformato in quello dellasomma di una serie. Il calcolo numerico dell’integrale si utilizza quando la funzione inte-granda e nota soltanto in un numero finito di punti o scaturisce da misure sperimentali,o ancora anche se e nota la primitiva, ma il suo calcolo puo essere cosı costoso che epreferibile utilizzare un metodo numerico per determinare una approssimazione. In tuttequeste situazioni e necessario approntare metodi numerici in grado di restituire un valoreapprossimato del calcolo dell’integrale.I metodi che vedremo consistono nell’approssimare o interpolare la funzione integrandamediante polinomi che risultano facilmente integrabili.

−3 −2 −1 0 1 2 30

5

10

15

20

25

30

35

40

45

50

Figure 2: f(x) = abs(x3 − 20) per x ∈ [−3, 3]

L’approssimazione del valore dell’integrale con la formula del punto medio, anche dettometodo dei rettangoli, porta ad individuare l’area del rettangolo in figura.

2011/01/12 26

−3 −2 −1 0 1 2 30

5

10

15

20

25

30

35

40

45

50

(a+b)/2 a b

Figure 3: IP M (f) = (b− a)f [(a + b)/2]

Formule di quadratura interpolatorie

Nella sua forma piu generale una formula di quadratura esprime un’approssimazionedi un integrale come una combinazione lineare dei valori della funzione integranda cioe:

∫ b

af(x)dx ≈

n∑

i=0

ωif(xi)

Sia Φ uno spazio di funzioni di dimensione n + 1 in cui sia sempre risolubile in modounico un problema di interpolazione su n + 1 punti distinti xj , j = 0 : n nell’intervallo[a, b] , allora e possibile trovare una base di funzioni per tale spazio, chiamate funzioniCARDINALI, per cui vale

ϕi(xj) =

{1 se i = j0 se i 6= j

Le funzioni cardinali risolvono il problema di interpolazione per i punti (xi, yi), i =0, 1, . . . , n, in modo banale, infatti la

Φ(x) =n∑

i=0

yiϕi(x)

verifica le condizioni di interpolazione.Osservazione

L’interpolante espressa nella base delle funzioni cardinali permette di esprimere banal-mente la soluzione di un problema di interpolazione analiticamente, ma non computazional-mente.Considerando l’interpolante nella base cardinale si esprime semplicemente una formula diquadratura interpolatoria, considerando yi = f(xi)

∫ b

af(x)dx '

∫ b

aΦ(x)dx

allora ∫ b

aΦ(x)dx =

∫ b

a

n∑

i=0

f(xi)ϕi(x)dx =

=n∑

i=0

f(xi)∫ b

aϕi(x)dx =

n∑

i=0

f(xi)ωi

2011/01/12 27

con

ωi =∫ b

aϕi(x)dx, i = 0, 1, . . . , n

Nel caso in cui consideriamo come spazio delle funzioni Φ lo spazio dei polinomi Pn,le funzioni cardinali sono i polinomi elementari di Lagrange

Li(x) = Πnj=0,j 6=i

(x− xj)(xi − xj)

Formule di quadratura di Newton-Cotes

Le formule di quadratura di Newton-Cotes si ottengono considerando n + 1 puntiequidistanti nell’intervallo chiuso [a, b] dati da:

xi = a + ih, i = 0, 1, . . . , n

con h = b−an per n > 0 e interpolando la funzione integranda mediante un polinomio

p ∈ Pn nei punti (xi, f(xi)), i = 0, 1, . . . , n. Se p(x) e il polinomio interpolante nella formadi Lagrange si avra

∫ b

af(x)dx '

∫ b

ap(x)dx =

n∑

i=1

f(xi)∫ b

aLi(x)dx

Quindi si procede all’integrazione dei polinomi Li(x), i = 0, 1, . . . , n per determinare icoefficienti ωi della formula di quadratura.Nell’integrare i polinomi elementari si esegue un cambiamento di variabile di integrazione,la nuova variabile t e tale che x = a + ht con t ∈ [0, n], poiche xi = a + ih, xj = a + jh edx = hdt avremo: ∫ b

aLi(x)dx

∫ b

aΠn

j=0,j 6=i

(x− xj)(xi − xj)

dx

=∫ n

0Πn

j=0,j 6=i

(a + ht− a− hj)(a + hi− a− hj)

hdt

= h

∫ n

0Πn

j=0,j 6=i

t− j

i− jdt

Osserviamo allora che i coefficienti ωi (anche detti pesi) sono dati da:

ωi =∫ n

0Πn

j=0,j 6=i

t− j

i− jdt

dipendono solo da n, non dipendono dalla funzione f che deve essere integrata e nemmenodagli estremi a e b dell’intervallo di integrazione.Le formule di Newton-Cotes piu comunemente usate sono quelle per

• n=1, Formula dei Trapezi

• n=2, Formula di Simpson

Ricaviamo queste formuleFormula dei Trapezi: n=1

ω0 =∫ 1

0

t− 10− 1

dt =∫ 1

0(1− t)dt = [t− t2

2] =

12

2011/01/12 28

ω1 =∫ 1

0

t− 01− 0

dt =∫ 1

0tdt = [

t2

2] =

12

da cui si ottiene ∫ b

af(x)dx '

∫ b

ap(x)dx =

h

2(f(a) + f(b))

La funzione f(x) e approssimata da una retta passante per i punti (x0, f(x0) e (x1, f(x1);cioe l’integrale e approssimativamete dato dall’area del Trapezio

−3 −2 −1 0 1 2 30

5

10

15

20

25

30

35

40

45

50abs(xf.3−20)

Formula di Simpson: n=2

ω0 =∫ 2

0

(t− 1)(0− 1)

(t− 2)(0− 2)

dt =12

∫ 2

0(t2 − 3t + 2)dt =

12[t3

3− t2

2+ 2t]20 =

12(83− 12

2+ 4) =

13

ω1 =∫ 2

0

(t− 0)(1− 0)

(t− 2)(1− 2)

dt = −∫ 2

0(t2 − 2t)dt =

= −(83− 4) =

43

ω2 =∫ 2

0

(t− 0)(2− 0)

(t− 1)(2− 1)

dt =12

∫ 2

0(t2 − t)dt =

12(83− 4

2) =

13

da cui ∫ b

af(x)dx '

∫ b

ap(x)dx =

h

3(f(x0) + 4f(x1) + f(x2))

La funzione f e approssimata da una parabola passante per (x0, f(x0), (x1, f(x1), (x2, f(x2)

EsempioApplichiamo la formula dei Trapezi e di Simpson al calcolo approssimato dell’integrale

di ∫ 1

0e−x2

dx

Il valore esatto alle prime 6 cifre significative e 0.746824.

• Trapezi∫ 1

0e−x2

dx ' h

2(e0 + e−1) =

12(1 + 0.367879) = 0.683939

2011/01/12 29

here

−3 −2 −1 0 1 2 30

5

10

15

20

25

30

35

40

45

50

• Simpsonin questo caso h = b−a

2 = 12

∫ 1

0e−x2

dx ' h

3(e0 + 4e−1/4 + e−1) =

16(1 + 4× 0.778801 + 0.367879) = 0.747180

−1.5 −1 −0.5 0 0.5 1 1.50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1y=exp(−x.2)

−1.5 −1 −0.5 0 0.5 1 1.50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1y=exp(−x.2) −− Simpson −−

Errore di Integrazione o grado di precisione di una formula di quadratura

Sia In una formula di quadratura che assumiamo essere esatta per ogni polinomio digrado ≤ q allora diciamo che In ha grado di precisone q. Il metodo dei trapezi e unaformula di quadratura con grado di precisione q = 1 cioe integra esattamente polinomi diprimo grado. ∫ b

af(x)dx =

h

2(f(x0) + f(x1)) + RT

∫ b

af(x)dx =

h

3(f(x0) + 4f(x1) + f(x2)) + RT

Dal teorema che fornisce l’errore di interpolazione polinomiale si ha che l’errore diintegrazione

RT =∫ b

a(x− x0)(x− x1) . . . (x− xn)

fn+1(ξ)(n + 1)!

dx

Nei casi particolari per n = 1 si ha:

• Errore nella formula dei Trapezi n=1

RT = − 112

h3f”(η) η ∈ [a, b]

2011/01/12 30

Dall’osservazione della formula dell’errore nel caso n = 1 si nota che il metodo dei Trapezie una formula di quadratura esatta per polinomi di grado ≤ 1, infatti la derivata seconda diquesti polinomi si annulla identicamente, e analogamnete la formula per n = 2 indica cheil metodo di Simpson sicuramente e esatta per polinomi di grado ≤ 2 cioe per polinomi digrado due perche hanno derivata terza nulla. Si dimostra pero che la formula di Simpsone esatta anche per polinomi di terzo grado cioe’ il grado di precisione e q ≤ 3.dimostrazione nel caso n=1

RT =∫ b

a(x− a)(x− b)

f”(ξ(x))(2)!

dx

ove ξ(x) ∈ [a, b]. Facciamo un cambiamento di variabile di integrazione, x = a + th; t ∈[0, 1]

RT = h

∫ 1

0(a + th− a)(a + th− a− h)

f”(ξ(a + th))(2)!

dt

=h3

2

∫ 1

0t(t− 1)f”(ξ(a + th))dt

applichiamo il teorena del valor medio per il calcolo dell’integrale:Teorema Se f(x) e una funzione continua, e g(x) e continua e non cambia segno nell’intervallo[a, b], allora esiste un punto ξ ∈ [a, b], tale che

∫ ba f(x)g(x)dx = f(ξ)

∫ ba g(x)dx

osserviamo che la parte t(t − 1) in [0, 1] non cambia segno, e se assumiamo che f” siacontinua, allora

RT =h3

2f”(ξ(a + th))

∫ 1

0t(t− 1)dt = −h3

12f”(η)

dove abbiamo posto η = ξ(a + th).

dimostrazione nel caso n=2 Per n=2 si ha la formula si Simpson che ha grado diprecisione q = 3 cioe integra esattamente polinomi di terzo grado. Infatti se consideriamoun polimnomio p3(x) lo possiamo scrivere nella forma

p3(x) = Ax3 + p2(x)∫ b

ap3(x)dx = A

∫ b

ax3dx +

∫ b

ap2(x)dx

sappiamo che per sua natura integra esattamente polinomi di secondo grado. Misuriamoallora l’errore di integrazione di questa formula di quadratura:

I(p3)− S2(p3) = A[I(x3)− S2(x3)] + [I(p2)− S2(p2)]

sappiamo che Simpson integra in modo esatto la parte quadratica del polinomio p3, quindila precisione risulta di grado 3 se integra in modo esatto x3, verifichiamolo:

S2(x3) =h

3[a3 + 4(

a + b

2)3 + b3] =

b− a

6[a3 +

12(a + b)3 + b3] =

=b− a

12[2a3 + (a + b)3 + 2b3] =

b− a

123[a3 + a2b + ab2 + b3] =

14[b4 − a4] = I(x3)

Quindi il metodo di quadratura di Simpson integra in modo esatto polinomi di terzo grado∫ ba p3(x)dx = S2(p3). Vale inoltre la seguente stima dell’errore di integrazione: se f ∈ CIV

[a,b]

allora esiste un punto η ∈ [a, b] tale che

RT = − 190

h5f IV (η) η ∈ [a, b]

2011/01/12 31

OsservazioneSe con s si indica il denominatore comune per i coefficienti frazionari ωi e indichiamo conσi = sωi, i = 0, 1, . . . , n i conseguenti numeri interi, allora la formula di Newton-Cotessi puo riscrivere come: ∫ b

af(x)dx ' b− a

ns

∑

i=0

σif(xi)

La seguente tabella riporta i coefficienti per le formule di newton-Cotes per n = 1, . . . , 8;si osservi che σi = σn−i per i = 0, . . . , n/2, per n ≥ 8 alcuni coefficienti assumono valorinegativi.

n σi ns Errore nome

1 1,1 2 −1/12h3f”(η) Trapezi2 1, 4, 1 6 −1/90h5f IV (η) Simpson3 1, 3, 3, 1 8 −3/80h5f IV (η) Cavalieri-Simpson4 7, 32, 12, 32, 7 90 −8/945h7fV I(η) Milne5 19, 75, 50, 50, 75, 19 288 . . . h7fV I(η)6 41 216 27 272. . . 840 . . . h9fV III(η) Weddle7 751 3577 1323 2989 . . . 17280 . . . h9fV III(η)8 989 5888 -928 10496 -4540 . . . 28350 . . . h11fX(η)

2011/01/12 32

here

−3 −2 −1 0 1 2 30

5

10

15

20

25

30

35

40

45

5073.0693

72.7344

trap

ezi c

ompo

sito

Formule composite Consistono nel suddividere l’intervallo [a, b] in sottointervalli ingenereale di uguale ampiezza [xi, xi+1] per i = 0, 1, . . . , m− 1, e in ciascuno di essi appli-care una formula di quadratura di grado basso.

Formula composita dei TrapeziSe si usa la formula dei Trapezi sul sottointervallo [xi, xi+1] si ottiene:

∫ xi+1

xi

f(x)dx =h

2[f(xi) + f(xi+1)]− h3

12f”(ηi)

con h = xi+1 − xi e xi < ηi < xi+1. L’errore di troncamento −h3

12f”(ηi) e detto errorelocale di trocamento

Per l’intervallo [a, b] suddiviso in m sottointervalli si ottiene

∫ b

af(x)dx =

m−1∑

i=0

∫ xi+1

xi

f(x)dx =

= h(12f(x0) + f(x1) + f(x2) + . . . + f(xm−1) +

12f(xm)) + RT

L’errore di troncamento globale RT e dato dalla somma degli errori locali cioe

RT = −h3

12

m−1∑

i=0

f”(ηi)

si noti che b− a = mh, e percio

RT = −b− a

12h2

∑m−1i=0 f”(ηi)

m

Se f” e continua in [a, b], allora esiste un η in questo intervallo, tale che

1m

m−1∑

i=0

f”(ηi) = f”(η)

da cui l’errore di troncamento totale per la formula dei Trapezi composita risulta

−b− a

12h2f”(η).

In particolare si osservi che aumentando m poiche h = b−am il valore di h diminuisce anzi

limh→0

RT = 0

2011/01/12 33

here

−3 −2 −1 0 1 2 30

5

10

15

20

25

30

35

40

45

5073.0693

Sim

pson

com

posi

to

72.7344

percio piu e fine la suddivisione dell’intervallo e migliore risulta l’approssimazione dell’integrale.

Formula composita di SimpsonSe m = 2k, cioe l’intervallo [a, b] e suddiviso in un numero pari di sottointervalli, alloral’integrale in ogni sottointervallo [x2i, x2i+2] puo essere calcolato con la formula di Simpson

∫ x2i+2

x2i

f(x)dx =h

3[f(x2i) + 4f(x2i+1) + f(x2i+2)]

−h5

90f IV (ηi)

con h = (x2i+2−x2i)2 e x2i < ηi < x2i+2. Per l’intero intervallo [a, b] si ottiene

∫ b

af(x)dx =

k−1∑

i=0

∫ x2i+2

x2i

f(x)dx =

h

3[f(x0) + 4f(x1) + 2f(x2) + . . .

+2f(x2i) + 4f(x2i+1) + 2f(x2i+2) + . . .+

+4f(x2k−1) + f(x2k) + RT

dove

RT = −h5

90

k−1∑

i=0

f IV (ηi)

Per le stesse argomentazioni usate per la formula dei Trapezi, l’errore si puo scrivere

RT = −b− a

180h4f IV (η)

dove η ∈ [a, b], sotto la condizione che f4 sia continua in [a, b].

EsempioSi calcola ∫ 1

0

11 + x

dx

usando le formule composite dei Trapezi e di Simpson. Scegliamo h = 1, .5, .25, 0.125 econfrontiamo i risultati con il valore esatto log 2 ' 0.693147 Si ottiene

2011/01/12 34

m h T(h) abs(T(h)-log(2)) S(h) abs(S(h)-log2)

1 1 0.750000 0.056853 - -2 0.5 0.708333 0.015186 0.694444 0.0012974 0.25 0.697024 0.003877 0.693254 0.0001078 0.125 0.694122 0.000975 0.693155 0.000008

Si noti che per la formula dei Trapezi poiche RT = O(h2), l’errore T(h)-log(2) dovrebbeessere diviso per un fattore 4 quando h si dimezza; analogamente per Simpson essendoRT = O(h4), l’errore S(h)-log2 dovrebbe essere diviso per 16 quando h si dimezza; i valoriin tabella rispecchiano questo andamento.

EsercizioSi determini il passo da utilizzare nella formula dei Trapezi composita, affinche la stimadell’integrale ∫ 1

0

1(x + 1)

dx

sia approssimata alla tolleranza 0.5× 10−2.

Formule Adattive

Il passo di integrazione h di una formula di quadrature puo essere scelto in modo dagarantire che l’errore sia inferiore ad una tolleranza ε > 0 prestabilita. Indichiamo con IA

il valore approssimato calcolato dell’integrale si vuole

|∫ b

af(x)dx− IA| ≤ ε

Occorre stabilire un criterio per determinare la stima dell’errore di integrazione, esistonodue approcci uno non adattativo e uno adattivo. In un approccio non adattivo i punti incui si valuta la funzione f(x) sono scelti senza tener conto del comportamento della fun-zione integranda. In questo caso se usassimo la formula di Simpson composita basterebberichiedere che

b− a

180h4maxx∈[a,b]|f IV (x)| < ε,

dove f4 indica la derivata quarta. Quindi basterebbe scegliere h sufficientemente piccolo inmodo da compensare l’apporto della derivata quarta. D’altra parte se f4 in valore assolutofosse grande solo in una piccola porzione dell’intervallo di integrazione si corre il rischiodi eseguire troppe valutazioni della f anche nelle parti in cui la funzione ha un compor-tamento lineare. In un approccio adattivo si utilizza una distribuzione non uniforme deisotto-intervalli nell’intervallo [a, b] in modo da garantire la stessa accuratezza della formulacomposita , ma con un numero inferiore di intervalli e quindi di valutazioni della f . A talfine serviranno uno stimatore dell’errore di quadratura ed una procedura che modifichi ilpasso di integrazione h conseguentemenete al soddisfacimento della tolleranza richiesta. Sisuddivide l’intervallo di integrazione in sottointervalli e si applica ricorsivamente a questiuna formula di quadratura, sfruttando una stima di arresto. Nella prossima sezione ve-dremo una tecnica nota come estrapolazione di Richardson che deriva da un risultato piugenerale da cui si puo estrarre una stima per l’errore di integrazione.

Estrapolazione di RichardsonE’ possibile determinare in modo automatico il numero m di sottointervalli in cui sud-dividere l’intero intervallo [a, b], per ottenere una stima accurata a meno di una certatolleranza. Il procedimento si basa sulla possibilita di stimare l’errore di integrazione RT .

2011/01/12 35

Cio puo essere fatto confrontando il valore dell’integrale ottenuto mediante la suddivisionein m sottointervalli, Im, con il valore I2m ottenuto con il doppio degli intervalli. Sia Iil valore esatto dell’integrale si puo dimostrare che per una formula di Newton-Cotes diprecisione s si ottiene la seguente stima dell’errore

|I − I2m| ' I2m − Im

2s+1 − 1

Si puo quindi innescare una procedura di raffinamenti successivi di [a, b] fino a che

I2m − Im

2s+1 − 1< tol

Caso Trapezi

Indichiamo con T (h) la formula dei trapezi composta con h = b−am ; allora vale

∫ b

af(x)dx− T (h) =

h2

12(f (1)(b)− f (1)(a))− h4

720(f (3)(b)− f (3)(a))+

h6

30240(f (5)(b)− f (5)(a)) + . . . + c2kh

2k(f2k−1(b)− f2k−1(a)) +O(h2k+2)

nell’ipotesi che la f sia derivabile in [a, b] almeno 2k + 2 volte.CorollarioAssumendo che la f(x) sia derivabile almeno 4 volte su [a, b], si ha:

∫ b

af(x)− T (

h

2) ' T (h/2)− T (h)

3+O(h4)

dimostrazione:dal teorema si ha: ∫ b

af(x)− T (h) = c2h

2 +O(h4)

e se il passo usato e h/2

∫ b

af(x)− T (h/2) = c2(h/2)2 +O(h4)

questo perche c2 e indipendente da h. Eliminando il termine in h2 fra le due equazioni, siottiene

4

(∫ b

af(x)− T (h/2)

)=

∫ b

af(x)dx− T (h) +O(h4)

da cui

3

(∫ b

af(x)− T (h/2)

)= T (h/2)− T (h) +O(h4)

come si voleva provare.Da questo si puo ricavare una formula che fornisce una approssima dell’integrale del quartoordine, cioe ∫ b

af(x)dx =

4T (h/2)− T (h)3

+O(h4)

2011/01/12 36

EsempioStimiamo ∫ 1

0e−x2

dx

usando la formula composita dei trapezi con h = 1, 0.5, 0.25. La stima viene poi miglioratamediante estrapolazione di Richardson. Il valore esatto alle prime 6 cifre significative e0.746824.

h T(h) Estrap. | ∫ ba fdx− T (h)| | ∫ b

a fdx− Estrap.|

1 0.683939 0.0628840.5 0.731370 0.747180 0.015453 0.000356420.25 0.742984 0.746855 0.003839 0.0000313797

Il corollario fornisce una stima del termine principale dell’errore. Si puo progettare unmetodo iterativo che dimezza il passo fino a che

|T (h/2)− T (h)3

| ≤ tol

quando questo si verifica sara anche∫ b

af(x)dx− T (h/2)| ≤ tol

Questo metodo iterativo puo essere applicato in modo adattivo.

2011/01/12 37

Caso SimpsonSe calcoliamo l’integrale con il metodo di Simpson con passo h = (b − a)/2 : S(h) e unSimpson composto con passo h/2:

∫ b

af(x)dx = S(h)− h5

90f (4)(η);

∫ b

af(x)dx = S(h/2)− b− a

180(h/2)4f (4)(η) =

= S(h/2)− 116

h5

90f (η).

procedendo

S(h)− S(h/2) ' h5

90f (4)(η)(1− 1

16)

da cui115

[S(h)− S(h/2)] ' 116

h5

90f (4)(η)

sostituendo nell’equazione precedente si ottiene

|∫ b

af(x)dx− S(h/2)| ' | 1

15[S(h/2)− S(h)] |

quindi c’e la possibilta di innescare un procedimento iterativo in cui si dimezza il passofino a che

| 115

[S(h/2)− S(h)] | ≤ tol

quando cio si verifica allora

|∫ b

af(x)dx− S(h/2)| ≤ tol

Metodo di Simpson AdattivoPer realizzare una routine che implementa un metodo adattivo si deve specificare l’intervallo[a, b], fornire una routine per la valutazione della f(x) per x ∈ [a, b] e scegliere una toller-anza tol. La routine deve calcolare un valore approssimato IA dell’integrale in modo che

|∫ b

af(x)dx− IA| ≤ tol

la routine deve anche essere in grado di determinare se la richiesta tolleranza non siaraggiungibile nei limiti definiti dal massimo numero di livelli di ricorsione. Durantel’esecuzione, ogni intervallo viene determinato per bisezione di un intervallo ottenutoprecedentemente durante il calcolo. Il numero effettivo di sottointervalli, cosı come laloro posizione e ampiezza dipende dalla f(x) e dalla tolleranza tol. Lo schema piu clas-sico applica ad ogni sottointervallo [xi, xi+1] la formula di Simpson semplice e la formuladi Simpson composta con m = 4, S(hi/2) ( che equivale a due Simpson semplici appli-cati rispettivamente a [xi, xi + hi/2] e [xi + hi/2, xi+1]. Entrambe S(hi) e S(hi/2) sonoapprossimazioni per

IAi '∫ xi+1

xi

f(x)dx.

Queste due approssimazioni servono per ottenere una stima dell’errore di integrazione.Se la tolleranza e raggiunta , S(hi/2) (o l’estrapolazione di Richardson relativa) vienepresa come valore dell’integrale su quell’intervallo. Se la tolleranza non viene raggiunta

2011/01/12 38

, il sottointervallo viene suddiviso a meta e il procedimento viene ripetuto su ognuno deisottointervalli piu piccoli.

Utilizzando il raffinamento automatico di Richardson, determinate il raffinamentodell’intervallo di definizione del seguente integrale utilizzando la formula dei Trapezi com-posita in modo da ottenere la seguente limitazione dell’errore

|RT | ≤ 1210−3

per il calcolo di

I =∫ 1

0e−x2

La function quad(fun, a, b) calcola una approssimazione dell’integrale della funzionefun utilizzando il metodo di Simpson adattivo: Dal manuale MatLab:

QUAD Numerically evaluate integral, adaptive Simpson quadrature.Q = QUAD(FUN,A,B) tries to approximate the integral of scalar-valuedfunction FUN from A to B to within an error of 1.e-6 using recursiveadaptive Simpson quadrature. FUN is a function handle. The functionY=FUN(X) should accept a vector argument X and return a vector resultY, the integrand evaluated at each element of X.

Q = QUAD(FUN,A,B,TOL) uses an absolute error tolerance of TOLinstead of the default, which is 1.e-6. Larger values of TOLresult in fewer function evaluations and faster computation,but less accurate results. The QUAD function in MATLAB 5.3 useda less reliable algorithm and a default tolerance of 1.e-3.

Q = QUAD(FUN,A,B,TOL,TRACE) with non-zero TRACE shows the valuesof [fcnt a b-a Q] during the recursion. Use [] as a placeholder toobtain the default value of TOL.

[Q,FCNT] = QUAD(...) returns the number of function evaluations.

Use array operators .*, ./ and .^ in the definition of FUNso that it can be evaluated with a vector argument.

Notes:Function QUADL may be more efficient with high accuracies and smoothintegrands.Function QUADV vectorizes QUAD for array-valued FUN.

Example:Q = quad(@myfun,0,2);

where myfun.m is the M-file function:%-------------------%function y = myfun(x)y = 1./(x.^3-2*x-5);%-------------------%

or, use a parameter for the constant:Q = quad(@(x)myfun2(x,5),0,2);

2011/01/12 39

where myfun2 is the M-file function:%----------------------%function y = myfun2(x,c)y = 1./(x.^3-2*x-c);%----------------------%

Class support for inputs A, B, and the output of FUN:float: double, single

See also quadv, quadl, dblquad, triplequad, trapz, function_handle.

Il metodo dei trapezi e’ invece implementato nella seguente function:

TRAPZ Trapezoidal numerical integration.Z = TRAPZ(Y) computes an approximation of the integral of Y viathe trapezoidal method (with unit spacing). To compute the integralfor spacing different from one, multiply Z by the spacing increment.

For vectors, TRAPZ(Y) is the integral of Y. For matrices, TRAPZ(Y)is a row vector with the integral over each column. For N-Darrays, TRAPZ(Y) works across the first non-singleton dimension.

Z = TRAPZ(X,Y) computes the integral of Y with respect to X usingthe trapezoidal method. X and Y must be vectors of the samelength, or X must be a column vector and Y an array whose firstnon-singleton dimension is length(X). TRAPZ operates along thisdimension.

Z = TRAPZ(X,Y,DIM) or TRAPZ(Y,DIM) integrates across dimension DIMof Y. The length of X must be the same as size(Y,DIM)).

Example: If Y = [0 1 23 4 5]

then trapz(Y,1) is [1.5 2.5 3.5] and trapz(Y,2) is [28];

Class support for inputs X, Y:float: double, single

See also sum, cumsum, cumtrapz, quad.

Lewis Fry Richardson (1881-1953) nacque a Newcastle in Gran Bretagna termino i suoistudi a King’s College a Cambridge nel 1903. Ricoprı varie cariche nell’industria e nell’universita,fu il primo a suggerire di usare tecniche matematiche per fare previsioni metereologiche ,mediante la risoluzione delle equazioni dei fluidi che dovrebbero governare la temperatura,la pressione atmosferica ecc. Fece questo negli anni venti, molto prima dello sviluppodei moderni computer ad alta velocita e fu a causa di questo che sviluppo la nozione dimetodi di estrapolazione per l’approssimazione numerica accurata di soluzioni basate suapprossimazioni meno precise.

2011/01/12 40

Esercizi di preparazione alla II prova parziale

1. Applica il metodo dei trapezi composto con passo h1 = 14 e h2 = 1

8 per fornire unastima del seguente integrale ∫ 1

0x(1− x2)dx

cerca di fornire una stima dell’errore di integrazione.

2. Dato il seguente integrale ∫ 1

−1

11 + x2

dx

prova a determinare il passo di integrazione che si deve usare per il metodo deitrapezi composto affinche l’errore di integrazione sia inferiore a 1

210−1

3. Poiche l’area del cerchio unitario e pari a π, segue che

π

2=

∫ 1

−1

√1− x2dx

Allora si puo approssimare π mediante una approssimazione di questo integrale.Applica il metodo di Simpson per calcolare valori approssimati di π commenta irisultati.