Interpolació polinòmica a partir de molts punts: splines cúbics
El problema de la interpolació
Sovint coneixem els valors d'una funció
en una taula de punts 
, o entorn d'ells, i a partir d'aquesta informació volem estimar el valor de en la resta de punts x del seu domini. En molts casos fins i tot podem tenir una funció desconeguda de la qual només podem mesurar algun valor.
en una taula de punts , o entorn d'ells, i a partir d'aquesta informació volem estimar el valor de en la resta de punts x del seu domini. En molts casos fins i tot podem tenir una funció desconeguda de la qual només podem mesurar algun valor. La interpolació consisiteix en aproximar aquesta funció per una funció polinòmica que passi pels punts que ja tenim a la taula. En la imatge que segueix, es veuen els punts que prenem (i en aquets cas la funció original en blau) i el polinomi interpolador (en magenta) i com seria l'error en un punt qualsevol. Es pot observar que pels punts originals, l'error és zero.
Fig1. Interpolació de 
(blau) per un polinomi de grau 4 (magenta). Error d'aproximació en x=0.63.
(blau) per un polinomi de grau 4 (magenta). Error d'aproximació en x=0.63.Una solució a aquest problema es presenta al problema 43 de la llista d'Àlgebra Lineal, però sobretot si hi ha molts punts i la funció no és polinòmica, pot ser una solució incorrecte.
no és polinòmica, pot ser una solució incorrecte. Una altra solució simple és traçar la linea poligonal en el pla pels punts
i aproximar per la funció 
que té aquesta poligonal com a gràfica. Però aquesta funció 
no és derivable en els vèrtexs de la poligonal, i aixó pot causar dificultats en el seu ús (canvis bruscs de velocitat que no són bons pels mecanismes).
per la funció que té aquesta poligonal com a gràfica. Però aquesta funció no és derivable en els vèrtexs de la poligonal, i aixó pot causar dificultats en el seu ús (canvis bruscs de velocitat que no són bons pels mecanismes).
Fig2. Interpolació poligonal de la funció 
en 8 punts. Falta de derivabilitat.
en 8 punts. Falta de derivabilitat.La solució
En aquesta pràctica demanarem com a punt de partida conéixer en una taula de punts 
els valors de i de la seva derivada 
:
els valors de i de la seva derivada : 
i aproximarem per un spline cúbic definit per aquesta taula. Un spline cubic serà doncs la funció 
definida a trossos en cada interval 
per polinomis 
respectivament, que compleixen les condicions:
per un spline cúbic definit per aquesta taula. Un spline cubic serà doncs la funció definida a trossos en cada interval per polinomis respectivament, que compleixen les condicions:- Cada polinomi
és de grau 3 o menys, 
(el polinomi te els valors de f en els extrems del seu interval),
(la derivada del polinomi te els valors de 
en els extrems del seu interval).
Si anomenem 
a cada interval entre dos punts on coneixem 
al polinomi de grau 3 que defineix el nostre spline, escrivint 
per a complir les condicions de spline cúbic explicades, els coeficients de 
en un interval qualsevol 
han de ser solució del sistema lineal:
a cada interval entre dos punts on coneixem al polinomi de grau 3 que defineix el nostre spline, escrivint per a complir les condicions de spline cúbic explicades, els coeficients de en un interval qualsevol han de ser solució del sistema lineal:
Trobar aquest polinomi s'anomena problema d'interpolació d'Hermite, i aquesta teoria matemàtica ens guaranteix que, si 
, el sistema lineal té solució única.
, el sistema lineal té solució única. Resolent el sistema anterior per a cada interval 
, obtenim els polinomis 
. La funció definida en cada interval pel seu polinomi respectiu és 
(te derivada 
contínua), perquè en els punts en que empalmem els polinomis de grau 3 ens hem assegurat de que tinguin el mateix valor i la mateixa derivada.
, obtenim els polinomis . La funció definida en cada interval pel seu polinomi respectiu és (te derivada contínua), perquè en els punts en que empalmem els polinomis de grau 3 ens hem assegurat de que tinguin el mateix valor i la mateixa derivada. Exemple:
Calculeu el polinomi d'Hermite que satisfà els valors de la taula:
que satisfà els valors de la taula:
Calculeu el valor interpolat per x=3.
Imposarem les condicions i resoldrem el sistema anterior
x0 = 0; xf = 5; %es poden posar dues instruccions a la mateixa línia
y0 = 1; yf = -2;
yp0 = 10; ypf = 2;
% matriu del sistema
A = [ x0^3 x0^2 x0 1;
3*x0^2 2*x0 1 0;
xf^3 xf^2 xf 1;
3*xf^2 2*xf 1 0];
% vector de termes independents
b = [y0; yp0; yf; ypf];
% solució del sistema
coef = A\b;
% els coeficients del polinomi: a3*x^3 + a2*x^2 + a1*x + a0
% son: a3 = coef(1); a2 = coef(2); a1 = coef(3); a0 = coef(4);
pHermite =[coef(1), coef(2), coef(3), coef(4)]; % es podria fer pHermite=coef directament
y3= polyval(pHermite,3) %aproximació de y(x) en el punt 3
y3 = 2.4160
% Noteu que la derivada s'aproximaria per la derivada del polinomi
derHermite = polyder(pHermite);
yp3 = polyval(derHermite,3)
yp3 = -4.3040
% verifiquem que en els punts de la taula, els valors coincideixen
y5= polyval(pHermite,5) %aproximació de y(x) en el punt 5
y5 = -2.0000
yp5 = polyval(derHermite,5)
yp5 = 2.0000
% Dibuixem el polinomi d'Hermite
xx=x0:0.1:xf; %ens calen molts punts per pintar la funció
yy=polyval(pHermite,xx);
plot(xx,yy)
grid
title('Polinomi d''Hermite')
La funció spline
Ara considerarem el mateix problema però per més d'un interval . Per a cada valor 
tindrem els valors de les seves 
respectivament. Per a cada interval, farem el càlcul dels coeficients dels polinomis d'Hermite respectius.
. Per a cada valor tindrem els valors de les seves respectivament. Per a cada interval, farem el càlcul dels coeficients dels polinomis d'Hermite respectius.El problema de calcular el valor interpolat en qualsevol altre punt dins l'interval 
es podria resoldre buscant el subinterval 
en el que està el punt i calcular el polinomi d'Hermite associat. Però si repetim el procés per molts punts llavors cada vegada que tinc un valor dins el mateix subinterval, hauria de refer el càlcul. En comptes d'això farem una vegada el càlcul del polinomi d'Hermite per cada subinterval i llavors aprofitarem els coeficients calculats.
es podria resoldre buscant el subinterval en el que està el punt i calcular el polinomi d'Hermite associat. Però si repetim el procés per molts punts llavors cada vegada que tinc un valor dins el mateix subinterval, hauria de refer el càlcul. En comptes d'això farem una vegada el càlcul del polinomi d'Hermite per cada subinterval i llavors aprofitarem els coeficients calculats.Al final d'aquest document teniu la funció spline_cubic que fa aquest càlcul.
Exemple:
Veiem un exemple de com muntar la funció spline i com trobar el valor interpolat en qualsevol nou punt.
% valors de la taula
x = [1,3,5,7]; % els intervals seran [1,3],[3,5],[5,7]
y = [0,4,9,12];
yp = [-1,2,5,10];
xn = 3.5; % nou punt a interpolar
% cridem a la funcio per a saber la taula de polinomis interpoladors del spline
taula_p = spline_cubic(x,y,yp);
% per saber en quin interval cau el punt es fa en dos passos:
% (1) busquem amb la funció find la primera component de x que es mes gran que xn
postk = find(x > xn, 1); % el primer valor és el 5 que està en el
% lloc (postk=3) del vector de les x
% (2) l'interval que conté xn es l'interval que
% acaba en el lloc postk (serà el segon interval)
k = postk-1;
% apliquem el polinomi interpolador k-essim per a estimar l'arrel cubica
p = taula_p(k,:);
yn = polyval(p,xn) %noteu que la solució yn=4.8750 és propera al valor en x=3
yn = 4.8750
Exercici:
Aplicarem aquesta funció a fer un programa que calculi aproximadament l'arrel cúbica real d'un nombre real positiu (la funció 
).
).Suposem que no sabem calcular arrels cúbiques però si sabem calcular potències enteres. Llavors, si construïm una taula a patir d'un vector de valors de y calculant les seves potències 
, la taula inversa [x, y] ens relaciona cada valor amb la seva arrel cúbica.
, la taula inversa [x, y] ens relaciona cada valor amb la seva arrel cúbica.Per exemple:
y = [1, 2, 3, 4, 5]
x = [1, 8, 27, 64, 125]
Mirant-ho al revés tenim que la 
.
.Per fer l'exercici, calcularem una taula de valors de la funció en 37 valors equiespaiats de y entre 1 i 10. D'aquesta manera tindrem informació sobre 37 valors de x entre 1 i 
pels que tindrem el valor de la seva arrel cúbica 
.
en 37 valors equiespaiats de y entre 1 i 10. D'aquesta manera tindrem informació sobre 37 valors de x entre 1 i pels que tindrem el valor de la seva arrel cúbica .La derivada 
en els punts de la taula ve donada pel teorema de la funció inversa, que us explicaran a Càlcul 1: com és la funció inversa de 
, la relació entre les derivades en un parell de valors 
qualsevol és que 
. Així obtenim el valorde la derivada 
pels 37 parells de la nostra taula.
en els punts de la taula ve donada pel teorema de la funció inversa, que us explicaran a Càlcul 1: com és la funció inversa de , la relació entre les derivades en un parell de valors qualsevol és que . Així obtenim el valorde la derivada pels 37 parells de la nostra taula.Ara passem aquesta taula de valors de 
a la nostra funcio spline_cubic.m per a que ens calculi la taula de polinomis que defineixen l'aproximació de per splines cúbics.
a la nostra funcio spline_cubic.m per a que ens calculi la taula de polinomis que defineixen l'aproximació de per splines cúbics. Un cop tenim la llista de polinomis taula_p que defineix el spline, ja podem usar-la per a estimar valors de 
: donat un x qualsevol en [1,1000] hem de localitzar entre quins 2 valors de la nostra taula de valors de x es troba, evaluar el polinomi de la fila corresponent de taula_p en x, i així obtenim el valor aproximat de .
: donat un x qualsevol en [1,1000] hem de localitzar entre quins 2 valors de la nostra taula de valors de x es troba, evaluar el polinomi de la fila corresponent de taula_p en x, i així obtenim el valor aproximat de . Comprovació:
Per a saber el marge d'error del nostre càlcul, podem comparar el valor de 
donat pel nostre spline amb el valor més precís (uns 14 decimals correctes) que dóna Matlab fent .
donat pel nostre spline amb el valor més precís (uns 14 decimals correctes) que dóna Matlab fent . Obs1: Només calcularem aquestes arrels cúbiques per a nombres reals 
. Amb això n'hi ha prou perquè tot nombre real es pot portar sobre aquest interval multiplicant-lo o dividint-lo per uns quants factors 1000, que contribueixen factors 10 a la seva arrel cúbica.
. Amb això n'hi ha prou perquè tot nombre real es pot portar sobre aquest interval multiplicant-lo o dividint-lo per uns quants factors 1000, que contribueixen factors 10 a la seva arrel cúbica.Obs2: El càlcul d'una funció a partir de taules de valors i interpolació per splines és recomanable en dispositius que tenen molta memòria per a emmagatzemar taules numèriques, però capacitat de bateria limitada per a fer càlculs. Si volem millorar la precissió del nostre càlcul, l'únic que s'ha de fer es refinar la taula de valors de x a partir de la que calculem el spline cúbic.
a partir de taules de valors i interpolació per splines és recomanable en dispositius que tenen molta memòria per a emmagatzemar taules numèriques, però capacitat de bateria limitada per a fer càlculs. Si volem millorar la precissió del nostre càlcul, l'únic que s'ha de fer es refinar la taula de valors de x a partir de la que calculem el spline cúbic.%------------------------------------------------------
% Aproximació dde l'arrel cúbica d'un nombre real
%------------------------------------------------------
% prenem 37 valors equiespaiats de y entre 1 i 10 (a distancia 0.25 entre ells)
y=linspace(1,10,37);
% valors de x son els de y al cub, calculats vectorialment
x=y.^3;
% Calculem y' segons el Tma de la funcio inversa, com s'explica
% a l'enunciat de la practica (y'=1/(3*y^2)), de manera vectorial
yp=1./(3*y.*y);
% cridem a la funcio per a saber la taula de polinomis interpoladors
% del spline agafant els vectors en l'ordre adequat
taula_p=spline_cubic(x,y,yp);
% interpolem en un punt qualsevol
xn=666; % volem l'arrel cubica d'aquest nombre
% Busquem amb find el lloc de la primera component
% de x que es mes gran que xn
postk=find(x>xn,1);
% l'interval que conte xn es l'interval que acaba en postk
k=postk-1;
% apliquem el polinomi interpolador k-essim per a estimar l'arrel cubica
p=taula_p(k,:);
yn=polyval(p,xn);
% error de la interpolacio: diferencia entre l'estimacio
% per interpolacio yn calculada i el valor amb 14 decimals
% correctes que calcula Matlab
% error absolut
error_abs=abs(yn-xn^(1/3))
error_abs = 8.4365e-08
% error relatiu
error_rel=error_abs/xn^(1/3)
error_rel = 9.6606e-09
% Obs1. Si hem reescalat el nombre xt per una potencia de 1000
% per a deixar-lo entre 1 i 1000, la mateixa potencia
% de 10 afectara al error absolut. En canvi l'error relatiu
% no varia al reescalar-lo.
%
% Obs2. Si volem tenir encara menys error: enlloc de 37 valors equiespaiats
% agafar-ne un nombre mes gran.
Codi de la funció spline_cubic:
Copieu aquest codi i guardeu-lo en un fitxer anomenat spline_cubic.m
function taula_p=spline_cubic(x,y,yp)
% Càlcul de la interpolació per splines cubics d'una funció
% a partir d'una taula de valors de la funció i de la seva derivada.
% INPUT: x=taula de punts on es sap el valor de la funció (ordenats
% de petit a gran, i en fila),
% y=taula de valors de la funcio en les abcises x (en fila),
% yp=taula de valors de la derivada de la funcio en les abcises x (en fila).
% OUTPUT: taula_p=matriu que conte en cada fila el polinomi interpolador
% de grau 3 del subinterval corresponent entre abcises.
%
% 2022/9/20
% comptem el nombre de valors a la taula (el nombre d'intervals serà un menys)
n = length(x);
taula_p = zeros(n-1,4); %taula per guardar els coeficients de cada polinomi d'Hermite
% bucle per cada interval entre valors
for k = 1:n-1
% extrems de l'interval on interpolarem
x0 = x(k);
xf = x(k+1);
y0 = y(k);
yf = y(k+1);
yp0 = yp(k);
ypf = yp(k+1);
% matriu de coeficients del sistema lineal
A = [x0^3 x0^2 x0 1;
3*x0^2 2*x0 1 0;
xf^3 xf^2 xf 1;
3*xf^2 2*xf 1 0];
% terme independent del sistema lineal
b = [y0; yp0; yf; ypf];
% trobem el polinomi interpolador d'Hermite (en columna)
p = A\b;
% i el guardem com a fila de la matriu
taula_p(k,:)=p';
end
end