Este pequeño programa muestra la diferencia entre interpolación y regresión

x=0.5:0.5:3;
y=[0.7,1.2,1.3,1.2,0.8,0.3];
hold on
plot(x,y,'ro','markersize',4,'markerfacecolor','r')

%regresión
p=polyfit(x,y,2); %polinomio de segundo grado
z=@(xx) polyval(p,xx);
fplot(z,[x(1),x(end)])

%interpolación
z=@(xx) spline(x,y,xx);
fplot(z,[x(1),x(end)])

hold off
grid on
legend('datos','regresión','interpolación')
xlabel('x')
ylabel('y');
title('Interpolación, regresión')

Interpolación lineal

La interpolación lineal es muy sencilla. Disponemos de pares de datos (x_k,y_k) k=1,2...n. Queremos conocer el valor de y para un valor cualesquiera de x en el intervalo x₁ a x_n. Supongamos que x está en el intervalo (x_k,x_k+1) tal como se muestra en la figura. Trazamos la recta que pasa por los puntos (x_k,y_k) y (x_k+1,y_k+1), cuya ecuación es

${\displaystyle y=\frac{(x-x_k)y_{k+1}-(x-x_{k+1})y_k}{(x_{k+1}-x_k)}}$

Creamos la función interpola_lineal para obtener el valor de y cuando se proporciona el valor de x.

function y0 = interpola_lineal(x,y,x0)
    y0=zeros(length(x0),1);
    for j=1:length(x0)
        indice=find(x>x0(j));
        k=indice(1)-1;
        y0(j)=((x0(j)-x(k))*y(k+1)-(x0(j)-x(k+1))*y(k))/(x(k+1)-x(k));
    end
end

Probamos la función para interpolar linealmente los siguientes pares de datos (0,0), (π/3, sin(π/3)), (2π/3, sin(2π/3)),(π, 0). Comprobamos que obtenemos los mismos resultado utilizando la función MATLAB interp1, pasándole la opción linear en el último parámetro.

x=0:pi/3:pi;  %datos
y=sin(x);
x0=[pi/6,pi/2,5*pi/6];  %interpolación lineal
y0=interpola_lineal(x,y,x0);
%y0=interp1(x,y,x0,'linear'); %esta es una función MATLAB
hold on
fplot('sin(x)',[0,pi])
plot(x,y)
plot(x,y,'bo','markersize',3,'markerfacecolor','b')
plot(x0,y0,'ro','markersize',4,'markerfacecolor','r')
set(gca,'XTick',0:pi/6:pi)
set(gca,'XTickLabel',{'0','\pi/6','\pi/3','\pi/2','2\pi/3','5\pi/6','\pi'})

hold off
grid on
xlabel('x');
ylabel('y')
title('Interpolación lineal')

Interpolación por un polinomio de grado n

Calculamos los coeficientes de un polinomio de tercer grado y=a₁x³+a₂x²+a₃x+a₄ que pase por los cuatro puntos. Obtenemos un sistema de cuatro ecuaciones con cuatro incógnitas.

${\displaystyle \begin{array}{l}\{\begin{array}{l}{a}_1{x}_1^3+a_2{x}_1^2+a_3{x}_1+a_4=y_1\\ a_1{x}_2^3+a_2{x}_2^2+a_3{x}_2+a_4=y_2\\ a_1{x}_3^3+a_2{x}_3^2+a_3{x}_3+a_4=y_3\\ a_1{x}_4^3+a_2{x}_4^2+a_3{x}_4+a_4=y_4\end{array}\\ \left(\begin{array}{cccc}{x}_1^3& x_1^2& x_1& 1\\ x_2^3& x_2^2& x_2& 1\\ x_3^3& x_3^2& x_3& 1\\ x_4^3& x_4^2& x_4& 1\end{array}\right)\left(\begin{array}{c}{a}_1\\ a_2\\ a_3\\ a_4\end{array}\right)=\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right)\end{array}}$

x=[0,pi/3,2*pi/3,pi];  %datos
y=sin(x);
xx=0:pi/90:pi;  %función
yy=sin(xx);
x0=[pi/6,pi/2,5*pi/6];  %interpolación lineal

A=vander(x);
p=A\y';  %coeficientes del polinomio

hold on
fplot('sin(x)',[0,pi])  % función
plot(x,y,'bo','markersize',3,'markerfacecolor','b')
z=@(xx) polyval(p,xx); %polinomio
fplot(z,[0,pi])
y0=polyval(p,x0); %valores interpolados
plot(x0,y0,'ro','markersize',4,'markerfacecolor','r')
set(gca,'XTick',0:pi/6:pi)
set(gca,'XTickLabel',{'0','\pi/6','\pi/3','\pi/2','2\pi/3','5\pi/6','\pi'})
hold off
xlabel('x');
ylabel('y')
grid on
title('Interpolación con un polinomio de tercer gardo')

Interpolación de Lagrange

Queremos encontrar los coeficientes de un polinomio de grado N

a₁x^N+a₂x^N-1+...+a_Nx+a_N+1

que pase por todos los pares de datos (x₁,y₁), (x₂,y₂), ...(x_N+1,y_N+1). Los coeficientes se pueden obtener resolviendo el sistema de ecuaciones:

${\displaystyle \{\begin{array}{l}{a}_1{x}_1^N+a_2{x}_1^{N-1}+a_N{x}_1+a_{N+1}=y_1\\ a_1{x}_2^N+a_2{x}_2^{N-1}+a_N{x}_2+a_{N+1}=y_2\\ \mathrm{.....}\\ a_1{x}_{N+1}^N+a_2{x}_{N+1}^{N-1}+a_N{x}_{N+1}+a_{N+1}=y_{N+1}\end{array}}$

En forma matricial

${\displaystyle \left(\begin{array}{ccccc}{x}_1^N& x_1^{N-1}& \mathrm{...}& x_1& 1\\ x_2^N& x_2^{N-1}& \mathrm{...}& x_2& 1\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ x_{N+1}^N& x_{N+1}^{N-1}& \mathrm{...}& x_{N+1}& 1\end{array}\right)\left(\begin{array}{c}{a}_1\\ a_2\\ \mathrm{...}\\ a_{N+1}\end{array}\right)=\left(\begin{array}{c}{y}_1\\ y_2\\ \mathrm{...}\\ y_{N+1}\end{array}\right)}$

Esta matriz se conoce con el nombre de Vandermonde y se construye con el comando vander(x), como podemos leer en el sistema de ayuda de MATLAB, A=vander(v) devuelve una matriz de Vandermonde cuyas columnas son las potencias del vector v, esto es A(i,j)=v(i)^(n-j)

Sea la siguiente tabla de datos tomada de una experiencia

Cuando el número N+1 de datos es relativamente pequeño, las primeras columnas de la matriz A pueden guardar números muy grandes, los efectos del redondeo pueden afectar al valor de los coeficientes a_i del polinomio.

Los valores del polinomio se puede obtener también utilizando la función polyfit que emplearemos más adelante en el ajuste por el procedimiento de mínimos cuadrados.

x=[0.97 1.12 2.92 3.00 3.33 3.97 6.10 8.39 8.56 9.44];
y=[2.58 0.43 0.06 5.74 7.44 8.07 6.37 2.51 1.44 0.52];
n=length(x); %número de pares de datos
%A=vander(x);
A=zeros(n); %líneas equivalentes a utilizar vander(x)
for j=1:n
    A(:,j)=(x').^(n-j);
end
p=A\y'; %sistema de ecuaciones lineales, y' es vector columna

%p=polyfit(x,y,n-1) %n-1 es el grado del polinomio
z=@(xx) polyval(p,xx);
hold on
plot(x,y,'bo','markersize',3,'markerfacecolor','b')
fplot(z,[x(1),x(end)])
xlabel('x')
ylabel('y')
grid on
title('Interpolación de Lagrange');
hold off

El script presenta opciones: La matriz A la obtenemos alternativamente, con el comando vander en vez de un bucle for. Los coeficientes del polinomio p, los obtenemos llamando a polyfit, en vez del operador división por la izquierda. Si llamamos a polyfit en la ventana de comandos aparece el siguiente mensaje

Warning: Polynomial is badly conditioned. Add points with distinct X
values, reduce the degree of the polynomial, or try centering and scaling
as described in HELP POLYFIT. 
> In polyfit (line 75)
  In xxx (line 11) 

p =   1.0e+04 *
0.0000 -0.0003 0.0057 -0.0595 0.3782 -1.4951 3.6430 -5.2142 3.9256 -1.1823

Un polinomio de Lagrange L_N(x) de grado N es la expresión

${\displaystyle \begin{array}{l}{L}_N(x)=\frac{(x-x_2)(x-x_3)\mathrm{...}(x-x_{N+1})}{(x_1-x_2)(x_1-x_3)\mathrm{...}(x_1-x_{N+1})}{y}_1+\\ \frac{(x-x_1)(x-x_3)\mathrm{...}(x-x_{N+1})}{(x_2-x_1)(x_2-x_3)\mathrm{...}(x_2-x_{N+1})}{y}_2+\mathrm{...}+\frac{(x-x_1)(x-x_2)\mathrm{...}(x-x_N)}{(x_{N+1}-x_1)(x_{N+1}-x_2)\mathrm{...}(x_{N+1}-x_N)}{y}_{N+1}\end{array}}$

Este polinomio pasa por todos los puntos (x_i,y_i).

Elaboramos una función denominada lagrange_3 que devuelva el valor interpolado yy de xx cuando se le pasa los vectores x e y que guardan las abscisas y ordenadas (x_i,y_i) de los datos

function yy = lagrange_3(x,y,xx)
   n = length(x); 
   for i=1:n 
       w(i)=prod(xx-x([1:i-1 i+1:n]))*y(i)/prod( x(i)-x([1:i-1 i+1:n]) ) ;
   end
   yy=sum(w);   
end

Llamamos a la función lagrange_3 para conocer los valores interpolados de xx=1.0 y xx=2.0

>> x=[0.97 1.12 2.92 3.00 3.33 3.97 6.10 8.39 8.56 9.44];
>> y=[2.58 0.43 0.06 5.74 7.44 8.07 6.37 2.51 1.44 0.52];
>> yy=lagrange_3(x,y,1.0)
yy =   10.0725
>> yy=lagrange_3(x,y,2.0)
yy = -203.7000

Alternativamente, elaboramos una función denominada lagrange_p que devuelva los coeficientes del polinomio. El polinomio de Lagrange es la suma de N+1 términos. El numerador de cada término es el producto de N binomios (x-x_i)

function p=lagrange_p(x,y)
    n=length(x); %n-1 es el grado del polinomio
    if length(y)~=n 
        error('x e y tienen que tener la misma longitud')
    end
    p=zeros(1,n);
    for i=1:n
        pol=[y(i)];
        for j=1:n
            if(i~=j)
%multiplica un polinomio por un binomio
                pol=conv([1 -x(j)],pol)/(x(i)-x(j)); 
            end
        end
        p=p+pol;
    end
end

Para obtener el producto de los binomios del numerador de cada término, utilizamos la función conv, véase la página Polinomios. Por ejemplo, queremos multiplicar los polinomios p₁=x³-2x-4 y p₂=x²+3x+4

>> p1=[1 0 -2 -4];
>> p2=[1 3 4];
>> p=conv(p1,p2)
p =     1     3     2   -10   -20   -16

Utilizamos la función lagrange_p para obtener los coeficientes del polinomio de Lagrange con los datos del ejemplo anterior.

>> x=[0.97 1.12 2.92 3.00 3.33 3.97 6.10 8.39 8.56 9.44];
>> y=[2.58 0.43 0.06 5.74 7.44 8.07 6.37 2.51 1.44 0.52];
>> p=lagrange_p(x,y)
p =  1.0e+004 *
0.0000 -0.0003 0.0057 -0.0595 0.3782 -1.4951 3.6430 -5.2142 3.9256 -1.1823

Obtenemos el mismo resultado, que resolviendo el sistema de N+1 ecuaciones lineales.

Una vez obtenido el polinomio de Lagrange calculamos el valor de y para valores de x que no están en la tabla.

>> xx=[1.0 2.0 3.5];
>> yy=polyval(p,xx)
yy =   10.0725 -203.7000    1.8338

Comparación

Como vamos a apreciar, la interpolación de Lagrange es un interesante ejercicio de programación, pero no poduce buenos resultados. Para probarlo, vamos a tomar nueve puntos de una función igualmente espaciados y vamos a comparar la gráfica de la función y la interpolación de los nueve puntos por el procedimiento de Lagrange. En la siguiente sección, explicamos un procedimiento mucho más exacto, los splines

Vamos a probar este procedimiento de interpolación, tomando nueve puntos equidistantes de la función

${\displaystyle y=\frac{\sin (x)}{x^2+1}}$

Ahora interpolamos mediante un polinomio de grado n-1 (ocho) que pasa por los n (nueve) puntos (x_i,y_i)

x=-4:4;
f=@(x) sin(x)/(1+x^2);
%dibuja los puntos
hold on
plot(x,y,'bo','markersize',3,'markeredgecolor','b','markerfacecolor','b')
n=length(x);
p=polyfit(x,y,n-1); %polinomio que pasa por los puntos
%p=lagrange_p(x,y); %se obtiene el mismo resultado
z=@(xx) polyval(p,xx);

fplot(z,[-4,4],'r')
fplot(f,[-4,4])
hold off
grid on
xlabel('x')
ylabel('y')
legend('función','interpolación')
title('Interpolación por Lagrange')

El resultado deja mucho que desear principalmente, en el primer y último intervalo.

Splines

Vamos a explicar mediante un ejemplo el procedimiento de interpolación denominado "splines", para generalizarlo después a cualquier conjunto de pares de datos.

Dado el conjunto de pares de datos (x₁,y₁), (x₂,y₂), (x₃,y₃), (x₄,y₄), (x₅,y₅) puntos de color rojo en la figura. Definimos la función S(x)

${\displaystyle \begin{array}{l}S(x)=\{\begin{array}{l}{s}_1(x)x_1\le x<x_2\\ s_2(x)x_2\le x<x_3\\ s_3(x)x_3\le x<x_4\\ s_4(x)x_4\le x<x_5\end{array}\\ s_i(x)=a_i{\left(x-x_i\right)}^3+b_i{\left(x-x_i\right)}^2+c_i\left(x-x_i\right)+d_{i}i=\mathrm{1,2,3,4}\end{array}}$

Cada una de las funciones s_i(x) en color azul en la figura, es un polinomio de tercer grado, cuya primera y segunda derivada es

${\displaystyle \begin{array}{l}\frac{d{s}_i}{dx}=3a_i{\left(x-x_i\right)}^2+2b_i\left(x-x_i\right)+c_i\\ \frac{d^2{s}_i}{d{x}^2}=6a_i\left(x-x_i\right)+2b_i\\ i=\mathrm{1,2,3,4}\end{array}}$

Para calcular las 4×4=16 incógnitas a₁,a₂,a₃,a₄, b₁,b₂,b₃,b₄, c₁,c₂,c₃,c₄, d₁,d₂,d₃,d₄, se imponen las siguientes condiciones:

1. La función S(x) pasa por todos los puntos (x_i,y_i)  i=1,2,3,4,5

${\displaystyle s_i(x_i)=y_i\{\begin{array}{l}{d}_1=y_1\\ d_2=y_2\\ d_3=y_3\\ d_4=y_4\\ a_4{\left(x_5-x_4\right)}^3+b_4{\left(x_5-x_4\right)}^2+c_4\left(x_5-x_4\right)+d_4=y_5\end{array}}$

2. La función S(x) es continua en x₂, x₃, x₄.

${\displaystyle s_{i-1}(x_i)=s_i(x_i)\{\begin{array}{l}{a}_1{\left(x_2-x_1\right)}^3+b_1{\left(x_2-x_1\right)}^2+c_1\left(x_2-x_1\right)+d_1=d_2\\ a_2{\left(x_3-x_2\right)}^3+b_2{\left(x_3-x_2\right)}^2+c_2\left(x_3-x_2\right)+d_2=d_3\\ a_3{\left(x_4-x_3\right)}^3+b_3{\left(x_4-x_3\right)}^2+c_3\left(x_4-x_3\right)+d_3=d_4\end{array}}$

3. La derivada primera de la función S(x) es continua en x₂, x₃, x₄.

${\displaystyle s{\text{'}}_{i-1}(x_i)=s{\text{'}}_i(x_i)\{\begin{array}{l}3a_1{\left(x_2-x_1\right)}^2+2b_1\left(x_2-x_1\right)+c_1=c_2\\ 3a_2{\left(x_3-x_2\right)}^2+2b_2\left(x_3-x_2\right)+c_2=c_3\\ 3a_3{\left(x_4-x_3\right)}^2+2b_3\left(x_4-x_3\right)+c_3=c_4\end{array}}$

4. Denominamos m₁, m₂, m₃, m₄, m₅ al valor de la derivada segunda de s_i(x) en los puntos x₁, x₂, x₃, x₄, x₅.

${\displaystyle s\text{'}{\text{'}}_i(x_i)=m_i\{\begin{array}{l}2b_1=m_1\\ 2b_2=m_2\\ 2b_3=m_3\\ 2b_4=m_4\\ 6a_4(x_5-x_4)+2b_4=m_5\end{array}}$

5. La derivada segunda de S(x) es continua en los puntos x₂, x₃, x₄

${\displaystyle s\text{'}{\text{'}}_{i-1}(x_i)=s\text{'}{\text{'}}_i(x_i)\{\begin{array}{l}6a_1(x_2-x_1)+2b_1=2b_2\\ 6a_2(x_3-x_2)+2b_2=2b_3\\ 6a_3(x_4-x_3)+2b_3=2b_4\end{array}}$

Supongamos que el espaciado entre dos puntos consecutivos es constante h=x_i-x_i-1, i=2,3,4,5.

Vamos a expresar las incógnitas a₁,a₂,a₃,a₄, b₁,b₂,b₃,b₄, c₁,c₂,c₃,c₄, d₁,d₂,d₃,d₄, en términos de h, el valor de las ordenadas y_i y el valor de la derivada segunda de S(x), m_i en cada punto x_i, i=1,2,3,4,5.

De (4) y (5) expresamos a₁,a₂, a₃, a₄ en términos de m₁, m₂, m₃, m₄, m₅.

${\displaystyle a_1=\frac{m_2-m_1}{6h}{a}_2=\frac{m_3-m_2}{6h}{a}_3=\frac{m_4-m_3}{6h}{a}_4=\frac{m_5-m_4}{6h}}$

Teniendo en cuenta (1) y (4), las ecuaciones (2) se expresan

${\displaystyle \{\begin{array}{l}\frac{m_2-m_1}{6h}{h}^3+\frac{m_1}{2}{h}^2+c_{1}h+y_1=y_2\\ \frac{m_3-m_2}{6h}{h}^3+\frac{m_2}{2}{h}^2+c_{2}h+y_2=y_3\\ \frac{m_4-m_3}{6h}{h}^3+\frac{m_3}{2}{h}^2+c_{3}h+y_3=y_4\end{array}}$

Que nos permite despejar c₁, c₂, c₃ en términos de m₁, m₂, m₃, m₄, y₁, y₂, y₃, y₄

${\displaystyle \begin{array}{l}{c}_1=\frac{y_2-y_1}{h}-\frac{m_2+2m_1}{6}h\\ c_2=\frac{y_3-y_2}{h}-\frac{m_3+2m_2}{6}h\\ c_3=\frac{y_4-y_3}{h}-\frac{m_4+2m_3}{6}h\end{array}}$

La última ecuación de (1) nos permite despejar c₄

${\displaystyle \begin{array}{l}\frac{m_5-m_4}{6h}{h}^3+\frac{m_4}{2}{h}^2+c_{4}h+y_4=y_5\\ c_4=\frac{y_5-y_4}{h}-\frac{m_5+2m_4}{6}h\end{array}}$

Las ecuaciones (3) se expresan

${\displaystyle \begin{array}{l}\{\begin{array}{l}3\frac{m_2-m_1}{6h}{h}^2+m_{1}h+\frac{y_2-y_1}{h}-\frac{m_2+2m_1}{6}h=\frac{y_3-y_2}{h}-\frac{m_3+2m_2}{6}h\\ 3\frac{m_3-m_2}{6h}{h}^2+m_{2}h+\frac{y_3-y_2}{h}-\frac{m_3+2m_2}{6}h=\frac{y_4-y_3}{h}-\frac{m_4+2m_3}{6}h\\ 3\frac{m_4-m_3}{6h}{h}^2+m_{3}h+\frac{y_4-y_3}{h}-\frac{m_4+2m_3}{6}h=\frac{y_5-y_4}{h}-\frac{m_5+2m_4}{6}h\end{array}\\ \{\begin{array}{l}{m}_1+4m_2+m_3=\frac{6}{h^2}\left(y_1-2y_2+y_3\right)\\ m_2+4m_3+m_4=\frac{6}{h^2}\left(y_2-2y_3+y_4\right)\\ m_3+4m_4+m_5=\frac{6}{h^2}\left(y_3-2y_4+y_5\right)\end{array}\end{array}}$

Tenemos tres ecuaciones y cinco incógnitas m₁, m₂, m₃, m₄, m₅

Fijamos los valores de la derivada segunda m₁ y m₅ en los puntos extremos. Supongamos que m₁=0 y m₅=0.

Hay otras posibilidades que se pueden consultar en el documento http://online.redwoods.edu/instruct/darnold/laproj/Fall98/SkyMeg/Proj.PDF

Despejamos m₂, m₃, m₄ del sistema de tres ecuaciones

${\displaystyle \begin{array}{l}\{\begin{array}{l}{m}_1+4m_2+m_3=\frac{6}{h^2}\left(y_1-2y_2+y_3\right)\\ m_2+4m_3+m_4=\frac{6}{h^2}\left(y_2-2y_3+y_4\right)\\ m_3+4m_4+m_5=\frac{6}{h^2}\left(y_3-2y_4+y_5\right)\end{array}\\ \left(\begin{array}{ccc}4& 1& 0\\ 1& 4& 1\\ 0& 1& 4\end{array}\right)\left(\begin{array}{c}{m}_2\\ m_3\\ m_4\end{array}\right)=\frac{6}{h^2}\left(\begin{array}{c}{y}_1-2y_2+y_3\\ y_2-2y_3+y_4\\ y_3-2y_4+y_5\end{array}\right)\end{array}}$

Caso general

En general, con n pares de datos tendremos el sistema

${\displaystyle \left(\begin{array}{cccccc}4& 1& 0& 0& \mathrm{...}& 0\\ 1& 4& 1& 0& \mathrm{...}& 0\\ 0& 1& 4& 1& \mathrm{...}& 0\\ 0& 0& 1& 4& \mathrm{...}& 0\\ \mathrm{....}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& 0& \mathrm{...}& 1& 4\end{array}\right)\left(\begin{array}{c}{m}_2\\ m_3\\ m_4\\ m_5\\ \mathrm{...}\\ m_{n-1}\end{array}\right)=\frac{6}{h^2}\left(\begin{array}{c}{y}_1-2y_2+y_3\\ y_2-2y_3+y_4\\ y_3-2y_4+y_5\\ y_4-2y_5+y_6\\ \mathrm{....}\\ y_{n-2}-2y_{n-1}+y_n\end{array}\right)}$

Fijamos los valores de los extremos del vector de las incógnitas m: m₁=0, m_n-1=0, a continuación, obtenemos mediante el operador división por la izquierda los valores de m₂, m₃, ...m_n-2, finalmente, se calculan:

• los coeficientes a₁,a₂,a₃,...a_n-1.

${\displaystyle a_1=\frac{m_2-m_1}{6h}{a}_2=\frac{m_3-m_2}{6h}{a}_3=\frac{m_4-m_3}{6h}\mathrm{....}{a}_{n-1}=\frac{m_n-m_{n-1}}{6h}}$

• los coeficientes b₁,b₂,b₃,...b_n-1.

${\displaystyle b_1=\frac{m_1}{2}{b}_2=\frac{m_2}{2}{b}_3=\frac{m_3}{2}\mathrm{....}{b}_{n-1}=\frac{m_{n-1}}{2}}$

• los coeficientes c₁,c₂,c₃,...c_n-1.

${\displaystyle \begin{array}{l}{c}_1=\frac{y_2-y_1}{h}-\frac{m_2+2m_1}{6}h\\ c_2=\frac{y_3-y_2}{h}-\frac{m_3+2m_2}{6}h\\ c_3=\frac{y_4-y_3}{h}-\frac{m_4+2m_3}{6}h\\ \mathrm{.....}\\ c_{n-1}=\frac{y_n-y_{n-1}}{h}-\frac{m_n+2m_{n-1}}{6}h\end{array}}$

• los coeficientes d₁,d₂,d₃,...d_n-1.

${\displaystyle d_1=y_1{d}_2=y_2{d}_3=y_3\mathrm{....}{d}_{n-1}=y_{n-1}}$

Codificación

En primer lugar, vamos a entender como trabaja la función diag de MATLAB. La función diag extrae vectores diagonales de la matriz A.

>> A=[1,2,3;4,5,6;7,8,9]
A =
     1     2     3
     4     5     6
     7     8     9
>> diag(A) %diagonal principal
ans =
     1
     5
     9
>> diag(A,1) %diagonal superior
ans =
     2
     6
>> diag(A,-1) %diagonal inferior
ans =
     4
     8

La función diag nos permite crear una matriz a partir de los vectores de sus diagonales

A=diag(4*ones(3,1))+diag(ones(2,1),1)+diag(ones(2,1),-1)
A =
     4     1     0
     1     4     1
     0     1     4

La función diff calcula la diferencia entre dos elementos consecutivos del vector y

${\displaystyle \begin{array}{l}y=[y_1,y_2,y_3,y_4,y_5]\\ d_1=diff(y)=[y_2-y_1,y_3-y_2,y_4-y_3,y_5-y_4]\\ d_2=diff(d_1)=[y_3-2y_2-y_1,y_4-2y_3-y_2,y_5-2y_4-y_3]\end{array}}$

El vector de los términos independientes se calcula aplicando dos veces la función diff al vector y de las ordenadas.

Estamos en condiciones de crear un script que calcule las incógnitas m₂, m₃, ... m_n-1, aplicando el operador división por la izquierda.

Vamos a probar este procedimiento de interpolación, tomando nueve puntos equidistantes de la función

${\displaystyle y=\frac{\sin (x)}{x^2+1}}$

x=-4:4;
f=@(x) sin(x)./(1+x.^2);
y=f(x);
h=x(2)-x(1); %espaciado constante

n=length(x); %número de pares de datos
%matriz de dimensión n-2
A=diag(4*ones(n-2,1))+diag(ones(n-3,1),1)+diag(ones(n-3,1),-1); 
s=diff(diff(y))*6/h^2; %vector de los términos indpendientes
mm=A\s'; %vector de las incógnitas
m=[0;mm;0]; %ceros en los extremos
a=zeros(1,n-1);
b=zeros(1,n-1);
c=zeros(1,n-1);
d=zeros(1,n-1);
for i=1:n-1
    a(i)=(m(i+1)-m(i))/(6*h);
    b(i)=m(i)/2;
    c(i)=(y(i+1)-y(i))/h-(m(i+1)+2*m(i))*h/6;
    d(i)=y(i);
end

hold on
plot(x,y,'bo','markersize',3,'markeredgecolor','b','markerfacecolor','b')
for i=1:n-1
    z=@(xx) a(i)*(xx-x(i)).^3+b(i)*(xx-x(i)).^2+c(i)*(xx-x(i))+d(i);
    fplot(z,[x(i),x(i+1)],'r')
end
fplot(f,[-4,4])
hold off
grid on
xlabel('x')
ylabel('y')
legend('función','interpolación')
title('Interpolación por splines')

con los nueve puntos unidos mediante polinomios de tercer grado vemos que la curva interpolada y la exacta están próximas.

Los puntos de color azul, son los nueve pares de datos, cada una de las curvas de color rojo entre dos puntos azules consecutivos es un polinomio de tercer grado que pasa por dichos puntos y cumple las condiciones de continuidad de su derivada primera y segunda.

Función spline de MATLAB

Utilizamos la función spline de MATLAB, pero es posible que las condiciones en los extremos que en el ejemplo anterior se han establecido en m₁=0, m_n-1=0, sean ahora diferentes.

x=-4:4;
f=@(x) sin(x)./(1+x.^2);
y=f(x);
%dibuja los puntos
hold on
plot(x,y,'bo','markersize',3,'markeredgecolor','b','markerfacecolor','b')

z=@(xx) spline(x,y,xx); %interpolación
fplot(z,[-4,4],'r')
fplot(f,[-4,4]) %función
hold off
grid on
xlabel('x')
ylabel('y')
title('Interpolación por splines de MATLAB')

Comparando ambas gráficas, vemos que el resultado que hemos obtenido creando nuestra propia función spline y utilizando la función spline de MATLAB es muy parecido

Aproximación de funciones mediante un polinomio

Muchas funciones matemáticas tales como sin(x), exp(x), se pueden representar cerca de x=0 por un desarrollo en serie de Taylor. Si el desarrollo en serie converge rápidamente, tomando los primeros términos se puede aproximar la función f(x) por un polinomio de grado n no muy elevado.

Puntos igualmente espaciados

Tomamos n valores de la función y_i=f(x_i) con i=1...n. y determinamos el polinomio de grado n-1 que pasa a través de los puntos (x_i,y_i). Sea la función

${\displaystyle y=\frac{1}{x^2+1}}$

Tomamos diez abscisas x_i espaciadas 0.8 en el intervalo [-4,4]. Representamos la función (en color rojo), los puntos (x_i,y_i) y el polinomio que pasa por dichos puntos (en color azul). Dado que la función es simétrica solamente representamos la parte positiva.

f=@(x) 1./(x.^2+1);
n=10;
a=-4;
b=4;

xx=linspace(a,b,n+1);
yy=f(xx);
p=polyfit(xx,yy,n);

x1=linspace(0,b,100);
y1=polyval(p,x1); %aproximación mediante el polinomio p
y2=f(x1);		%función
hold on
plot(x1,y1,'b')
plot(x1,y2,'r')
plot(xx(xx>=0),yy(xx>=0),'ro','markersize',3,'markerfacecolor','r')
hold off
xlabel('x')
ylabel('y')
title('Aproximación de una función')

Puntos espaciados mediante la fómula de Chebyshev

Sorprendentemente, si los puntos no están igualmente espaciados, obtenemos una mejor aproximación. Si tomamos los puntos en el intervalo [a,b] espaciados de acuerdo a la fórmula denominada puntos de Chebyshev.

${\displaystyle x_i=\frac{a+b}{2}+\frac{a-b}{2}\cos \left(\frac{\pi }{n}\left(i-\frac{1}{2}\right)\right)1\le i\le n}$

La función linspace, nos crea un vector de abscisas x_i igualmente espaciados entre a y b. La función que hemos denominado lincheby crea un vector de abscisas x_i, espaciados de acuerdo a la fórmula anterior.

f=@(x) 1./(x.^2+1);
lincheby=@(a,b,n) (a+b)/2+(a-b)/2*cos(pi/n*(1/2:n));
n=10; 
a=-4;
b=4;

xx=lincheby(a,b,n+1);
yy=f(xx);
p=polyfit(xx,yy,n); 

x1=linspace(0,b,100);
y1=polyval(p,x1); %aproximación mediante el polinomio p
y2=f(x1);		  %función
hold on
plot(x1,y1,'b')
plot(x1,y2,'r')
plot(xx(xx>=0),yy(xx>=0),'ro','markersize',3,'markerfacecolor','r')
hold off
xlabel('x')
ylabel('y')
title('Aproximación de una función')

Como podemos apreciar la mejora es importante