Laboratorio 4
Interpolazione polinomiale
Esercizio 1
Per agevolare la risoluzione dell’esercizio digitiamo tutti i comandi necessari in un file che
salveremo con estensione .sce (si ricorda che per far eseguire i comandi da Scilab si digita
da linea di comando exec nomefile.sce). Per analizzare il comportamento del polinomio
interpolatore con diverso grado, dovremo volta per volta aggiornare soltanto la variabile che
contiene il grado del polinomio.
1. Iniziamo realizzando il grafico della funzione f (x):
f=’sin(x).*x’;
a=-1;
b=3;
x_disegno=[a:0.01:b];
x=x_disegno;
scf(1);
plot2d(x_disegno,eval(f));
2. Per calcolare il polinomio interpolatore di Lagrange utilizziamo la funzione polyfit.sci.
La funzione polyfit.sci richiede come parametri in ingresso: il vettore contenente i
nodi di interpolazione, il vettore contenente i valori della funzione da interpolare su tali
nodi ed infine il grado del polinomio interpolatore da calcolare (che nel caso in esame
sarà pari al numero di nodi meno 1 in quanto si vuole un’interpolazione di Lagrange).
In uscita la funzione ci ritorna il polinomio interpolatore espresso in forma simbolica
(del tipo a0 + a1 ∗ X + a2 ∗ X 2 + ...... + an+1 ∗ X n ).
Definiamo i nodi su cui verrà costruito il polinomio interpolatore e valutiamo la funzione
f su tali nodi:
n=2;
h=(b-a)/n;
xnodi=[a:h:b];// vettore contenente n+1 nodi equispaziati
x=xnodi;
1
ynodi=eval(f);
grado=(max(size(xnodi)))-1; //grado del polinomio interpolatore
Carichiamo la funzione polyfit.sci e calcoliamo il polinomio interpolatore:
getf polyfit.sci
P=polyfit(xnodi,ynodi,grado);
Utilizziamo la funzione horner che valuta il polinomio interpolatore su un vettore di
punti passato in ingresso alla funzione.
La funzione richiede come primo parametro d’ingresso il polinomio interpolatore in
forma simbolica (tale variabile é il parametro d’uscita della funzione polyfit) e come
secondo parametro d’ingresso il vettore di punti su cui vogliamo valutare il polinomio
interpolatore.
y_pol=horner(P,x_disegno);
Alla funzione horner abbiamo passato P, parametro che ci restituisce polyfit.sci,
come primo parametro; come secondo parametro gli passiamo lo stesso insieme di punti
utilizzato per visualizzare la funzione f .
Rappresentiamo graficamente la funzione, il polinomio interpolatore calcolato e i valori
assunti dalla funzione f nei nodi d’interpolazione:
scf(2)
plot2d([x_disegno’ x_disegno’],[y_disegno’ y_pol’],style=[2 5]);
plot2d(xnodi,ynodi,style=-4);
In figura 1 riportiamo la funzione e i polinomi calcolati. Si noti che per n=5 non si
riesce a distinguere visivamente la funzione e il polinomio interpolatore.
3. Visualizziamo l’errore ε e calcoliamo l’errore massimo commesso:
errore=abs(y_disegno-y_pol);
scf(3);
plot2d(x_disegno,errore,style=5);
e_max=max(errore);
In figura 2 vediamo l’andamento dell’errore mentre εmax vale:
per n=2
2
2.2
♦
1.8
♦
♦
1.4
1.0
♦
♦
0.6
♦
♦
0.2
♦
♦
-0.2
-1.0
-0.6
-0.2
0.2
0.6
1.0
1.4
1.8
2.2
2.6
3.0
Figura 1: Funzione e polinomi interpolanti per n=2 (in rosso), per n=3 (in azzurro), per n=5
(in verde)
--> e_max
=
1.14957
per n=3
--> e_max
=
0.26272
per n=5
--> e_max
=
0.0201661
I valori numerici assunti dall’errore massimo per i tre valori di n ci confermano che
l’interpolazione risulta sempre più accurata all’aumentare del grado del polinomio.
3
1.2
1.0
0.8
0.6
0.4
0.2
0.0
-1.0
-0.6
-0.2
0.2
0.6
1.0
1.4
1.8
2.2
2.6
3.0
Figura 2: Andamento dell’errore nel caso n=2 (in rosso), n=3 (in azzurro), n=5 (in verde)
Esercizio 2
1. Costruiamo l’interpolante Πn f di grado n=5,10,20 della funzione f considerando nodi
equispaziati sull’intervallo I. Per ciascun valore di n si vuole rappresentare graficamente
Πn f e l’errore |f (x) − Π(x)|.
Come per l’esercizio precedente digitiamo i comandi in un file che salveremo con estensione .sce.
In primo luogo definiamo la funzione di Runge:
f_Runge=’(1)./(1+x.^2)’;
Definiamo gli estremi dell’intervallo I nelle variabili a,b e il vettore contenente n+1 nodi
equispaziati nell’intervallo di interesse (xnodi).
n=5;
a=-5;
b=5;
h=(b-a)/n; //distanza tra i nodi
xnodi=[a:h:b]; //vettore contenente n+1 nodi equispaziati
4
Procediamo valutando la funzione in corrispondenza dei nodi prima definiti:
x=xnodi;
ynodi=eval(f_Runge);
Per realizzare il polinomio interpolatore utilizzamo le due funzioni polyfit.sci e
horner in analogia a quanto fatto nell’esercizio 1.
grado=(max(size(xnodi)))-1; //grado del polinomio interpolatore
getf polyfit.sci //caricamento della funzione
P=polyfit(xnodi,ynodi,grado);
x_disegno=[a:0.01:b];
y_pol=horner(P,x_disegno);
Realizziamo il grafico della funzione nell’intervallo I e del polinomio interpolatore Π n f .
x=x_disegno;
y_disegno=eval(f_Runge);
scf(1)
plot2d([x_disegno’ ,x_disegno’],[y_disegno’, y_pol’],style=[2 5]);
plot2d(xnodi,ynodi,style=-4);
Infine analizziamo l’andamento dell’errore definito come |f (x) − Πn (x)| e rappresentiamolo graficamente:
scf(2)
errore=abs(y_disegno-y_pol);
plot2d(x_disegno,errore,style=5);
I risultati sono mostrati nelle figure 3, 4 e come si vede l’errore di interpolazione risulta
maggiore ai bordi dell’intervallo rispetto ai valori assunti al centro. Tale comportamento
si accentua sempre di piú all’aumentare di n (numero dei nodi equispaziati).
Si verifica che esistono dei punti x interni all’intervallo tali per cui:
lim |f (x) − Πn f (x)| 6= 0
n→∞
non si ha convergenza uniforme dell’errore.
5
2. Consideriamo ora i nodi di Chebyshev, ed usiamo le seguenti istruzioni in Scilab per
costruire il polinomio di interpolazione sui nodi x̂i :
k=[0:n];
t=-cos(%pi*k/n);//vettore contenete i nodi di interpolazione
xnodi=((b-a)/2)*t+(a+b)/2; //riscalatura
Calcoliamo il polinomio interpolatore e rappresentiamo la funzione e l’interpolante come
nel caso dei nodi equispaziati:
x=xnodi;
ynodi=eval(f_Runge);
grado=(max(size(xnodi)))-1;
getf polyfit.sci
P=polyfit(xnodi,ynodi,grado);
y_pol=horner(P,x_disegno);
scf(3);
plot2d([x_disegno’ x_disegno’],[y_disegno’ y_pol’],style=[2 5]);
plot2d(xnodi,ynodi,style=-4);
Valutiamo anche l’errore di interpolazione e rappresentiamo l’andamento su un grafico:
scf(4)
errore=abs(y_disegno-y_pol);
plot2d(x_disegno,errore,style=5);
I risultati riportati nelle figure 5,6 mostrano come l’interpolazione di Chebyshev garantisca una distribuzione uniforme dell’errore sull’intervallo; precisamente si osserva che,
per n = 10, l’errore massimo è di ordine 0.1 mentre nel caso dei nodi equispaziati era
circa pari a 2. L’accuratezza dell’interpolazione aumenta al crescere di n come si può
facilmente verificare.
6
2.1
1.7
1.3
♦
0.9
♦
0.5
♦
♦
0.1 ♦
♦
♦
♦
♦
♦
♦
4
5
-0.3
-5
-4
-3
-2
-1
0
1
2
3
Figura 3: Funzione ed interpolante nel caso n=5 (in rosso) e n=10 (in azzurro) su nodi
equispaziati
2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
-5
-4
-3
-2
-1
0
1
2
3
4
5
Figura 4: Errore nel caso n=5 (in rosso) e n=10 (in azzurro) su nodi equispaziati
7
♦
1.0
0.9
0.8
0.7
0.6
0.5
0.4
♦
0.3
♦
0.2
0.1
♦
♦♦
♦
-5
-4
♦
♦
♦♦
4
5
0.0
-3
-2
-1
0
1
2
3
Figura 5: Funzione ed interpolante nel caso n=5 (in rosso) e n=10 (in azzurro) su nodi di
Chebyshev
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
-5
-4
-3
-2
-1
0
1
2
3
4
5
Figura 6: Errore nel caso n=5 (in rosso) e n=10 (in azzurro) su nodi di Chebyshev
8