% Activate the following line by filling in the right side. If for example the name of the root file is Main.tex, write
% "...root = Main.tex" if the chapter file is in the same directory, and "...root = ../Main.tex" if the chapter is in a subdirectory.
 
%!TEX root = ../Metodos_Numericos.tex  

\chapter[Ajuste de curvas]{Ajuste de curvas}

\section{Regresión lineal por el método de mínimos cuadrados}

El ejemplo más simple de una aproximación por mínimos cuadrados es mediante el ajuste de un conjunto de pares de observaciones: $[x_1, y_1]$, $[x_2, y_2]$, $[x_3, y_3]$,..., $[x_n, y_n]$ a una línea recta. La expresión matemática de la línea recta es:

$$
y = a_0 + a_1x + e
$$
 
donde  $a_0$  y  $a_1$, son los coeficientes que representan el cruce con el eje y y la pendiente de la línea y e representa el error de nuestra aproximación. En la figura \ref{fig:pros}, se muestra un ejemplo de puntos a los que se les desea ajustar a una linea recta.

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_6/figuras/puntos.eps, width=10.0cm}}
   \caption{Un conjunto de puntos}
   \label{fig:ptos}
\end{figure}

Una estrategia, para ajustar a la mejor línea, es minimizar la suma al cuadrado de los errores para todos los datos disponibles

\begin{equation*}
E(a) = \sum_{i=1}^{n} e_i^2 = \sum_{i=1}^{n} (y_i - a_0 - a_1 x_i)^2
\end{equation*}


Esta misma ecuación la podemos escribir en forma matricial como

\begin{equation*}
E(a) = [e_1, e_2, \cdots, e_n]
\left[
\begin{array}{c}
  e_1   \\
  e_2   \\
  \vdots\\
  e_n   
\end{array}
\right]
\end{equation*}

Los elementos del vector de error $e_i$ en forma matricial queda como

\begin{equation*}
e_i = y_i - [1, x_i]
\left[
\begin{array}{c}
  a_0  \\
  a_1     
\end{array}
\right]
\end{equation*}


y en general

\begin{equation*}
\left[
\begin{array}{c}
  e_1   \\
  e_2   \\
  e_3 \\
  \vdots\\
  e_n   
\end{array}
\right] =
\left[
\begin{array}{c}
  y_1   \\
  y_2   \\
  y_3 \\
  \vdots\\
  y_n   
\end{array}
\right] -
\left[
\begin{array}{cc}
  1 & x_1   \\
  1 & x_2   \\
  1 & x_3 \\
  \vdots & \vdots\\
  1 & x_n   
\end{array}
\right]
\left[
\begin{array}{c}
  a_0   \\
  a_1   
\end{array}
\right]
\end{equation*}

En forma compacta tenemos que $e = y – Ma$, donde $M$ es la matriz de coordenadas en $x$ y $a$ el vector de parámetros.

\subsection{Ajuste por mínimos cuadrados}

Si queremos encontrar el vector de parámetros a que minimiza nuestra suma de cuadrados, tenemos que calcular la derivada de la función de error respecto al vector de parámetros e igualar a cero.

\begin{eqnarray*}
\frac{\partial E(a)}{\partial a} = \frac{\partial}{\partial a} [y - Ma]^T [y-Ma] =0\\
\frac{\partial E(a)}{\partial a} = -2M^T [y - Ma] =0
\end{eqnarray*}

El valor del vector de parámetros a los calculamos resolviendo el siguiente sistema de ecuaciones

\begin{equation*}
[M^TM]a = M^Ty
\end{equation*}

Es común encontrar la solución de este sistema como:

\[
\left[
\begin{array}{c|c}
  \sum_{i=1}^{n}1    &  \sum_{i=1}^{n}x_i  \\
  \hline
  \sum_{i=1}^{n}x_i &   \sum_{i=1}^{n}x_i^2   
\end{array}
\right]
\left[
\begin{array}{c}
  a_0  \\
  \hline
  a_1   
\end{array}
\right] =
\left[
\begin{array}{c}
  \sum_{i=1}^{n}y_i  \\
  \hline
  \sum_{i=1}^{n}x_i y_i 
\end{array}
\right]
\]

La implementación en Matlab es:

\begin{verbatim}
function a = AjustePolinomial (x, y, grado) 

    N = length(x);
    M = ones(N, 1);

    for k=1:grado
        M = [M,x.^k];
    end

    a=inv(M'*M)*M'*y;
end
\end{verbatim}

 
\subsection{Ejemplo 1}

Hacer el ajuste a una línea recta de los siguientes valores

\begin{center}
\begin{tabular}{r|r}
$x$ & $y$\\
\hline
1.00 & 0.50\\
2.00 & 2.50\\
3.00 & 2.00\\
4.00 & 4.00\\
5.00 & 3.50\\
6.00 & 6.00\\
7.00 & 5.50\\
\end{tabular}
\end{center}
 
La matriz $M$ queda como

\[
M = 
\left[
\begin{array}{cc}
1 & 1.00\\
1 & 2.00\\
1 & 3.00\\
1 & 4.00\\
1 & 5.00\\
1 & 6.00\\
1 & 7.00\\
\end{array}
\right]
\]

\[
M^T M= 
\left[
\begin{array}{ccccccc}
1 & 1 & 1 & 1 & 1 & 1 & 1\\ 
1.00 & 2.00 & 3.00 & 4.00 & 5.00 & 6.00 & 7.00\\
\end{array}
\right]
\left[
\begin{array}{cc}
1 & 1.00\\
1 & 2.00\\
1 & 3.00\\
1 & 4.00\\
1 & 5.00\\
1 & 6.00\\
1 & 7.00\\
\end{array}
\right] =
\left[
\begin{array}{rr}
  7.00 & 28.00  \\
28.00 & 140.00\\ 
\end{array}
\right]
\]

\[
M^T y= 
\left[
\begin{array}{ccccccc}
1 & 1 & 1 & 1 & 1 & 1 & 1\\ 
1.00 & 2.00 & 3.00 & 4.00 & 5.00 & 6.00 & 7.00\\
\end{array}
\right]
\left[
\begin{array}{cc}
0.50\\
2.50\\
2.00\\
4.00\\
3.50\\
6.00\\
5.50\\
\end{array}
\right] = 
\left[
\begin{array}{r}
24.00 \\
119.50
\end{array}
\right]
\]

Aplicando las formulas anteriores, tenemos que el sistema de ecuaciones a resolver es

 \[
\left[
\begin{array}{rr}
  7.00 & 28.00  \\
28.00 & 140.00\\ 
\end{array}
\right]
\left[
\begin{array}{r}
a_0 \\
a_1\\
\end{array}
\right] = 
\left[
\begin{array}{r}
24.00 \\
119.50
\end{array}
\right]
\]
La solución es $a = [0.07142, 0.8392]$ y en la figura  \ref{fig:ptos_linea} se muestra el ajuste encontrado

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_6/figuras/puntos_linea.eps, width=10.0cm}}
   \caption{Un conjunto de puntos y la linea recta ajustada}
   \label{fig:ptos_linea}
\end{figure}

Para ejecutar dar
\begin{verbatim}
AjustePolinomial([1,2,3,4,5,6,7]', [0.5, 2.5, 2, 4, 3.5, 6, 5.5]',1)
\end{verbatim}

\subsection{Regresión polinomial}

Podemos generalizar el caso de la regresión lineal y extenderla a cualquier polinomio de orden $m$, hacemos la siguiente representación.
\begin{equation*}  
\left[
\begin{array}{ccccccc}
  1 & x_1  & x_1^2  & x_1^3 & \cdots & x_1^m \\
  1 & x_2  & x_2^2  & x_2^3 & \cdots & x_2^m \\
  1 & x_3 & x_3^2  & x_3^3& \cdots & x_3^m \\
  \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\
  1 & x_n  & x_n^2  & x_n^3 & \cdots & x_n^m 
\end{array}
\right]
\left[
\begin{array}{c}
a_0\\
a_1\\
a_2\\
a_3\\
\vdots\\
a_m
\end{array}
\right]
=
\left[
\begin{array}{c}
y_0\\
y_1\\
y_2\\
\vdots\\
y_n
\end{array}
\right]
\end{equation*}

La cual podemos escribir de manera compacta como

\begin{equation*}
M a = y
\end{equation*}

Si pre-multiplicamos ambos lados de la ecuación tenemos

\begin{equation}
[M^TM]a = M^Ty
\end{equation}

El sistema se resuelve para $a$ o simplemente se calcula haciendo
 
\begin{equation*}
a = [M^TM]^{-1}[M^Ty]
\end{equation*}

\subsection{Ejemplo 2}

Ajustar a un polinomio de segundo orden los datos en la siguiente tabla.

\begin{center}
\begin{tabular}{r|r}
$x$ & $y$ \\
\hline
0.00 & 2.10\\
1.00 & 7.70\\
2.00 & 13.60\\
3.00 & 27.20\\
4.00 & 40.90\\
5.00 & 61.10\\
\end{tabular}
\end{center}

En este caso por ser un polinomio de segundo grado tenemos una matriz $M$ dada por
\[
M = 
\left[
\begin{array}{ccc}
     1  &   0  &   0\\
     1  &   1  &   1\\
     1  &   2  &   4\\
     1  &   3  &   9\\
     1  &   4  &  16\\
     1  &   5  &  25\\
\end{array}
\right]
\]

El sistema de ecuaciones que debemos resolver es:

\begin{equation*}
\left[
\begin{array}{rrr}
  6.00 &   15.00   &  55.00\\
15.00 &   55.00  & 225.00\\
55.00 & 225.00  & 979.00
\end{array}
\right]
\left[
\begin{array}{r}
a_0\\
a_1\\
a_2
\end{array}
\right]
=
\left[
\begin{array}{r}
      152.60\\
      585.60\\
    2488.80
\end{array}
\right]
\end{equation*}
 
La solución del sistema es $a = [2.4785, 2.3592, 1.8607]$ y el ajuste da como resultado el polinomio $p(x) = 2.4785 + 2.3592x +  1.8607x2$. En la figura \ref{fig:ptos_2g} se muestra la aproximación a un polinomio de segundo orden. 

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_6/figuras/puntos_2g.eps, width=10.0cm}}
   \caption{Un conjunto de puntos y su aproximación cuadrática}
   \label{fig:ptos_2g}
\end{figure}

Para correr en Matlab hacer
\begin{verbatim}
AjustePolinomial([0,1,2,3,4,5]', [2.1, 7.70, 13.60, 27.20, 40.90, 61.10]',2)
\end{verbatim}

\section{Interpolación lineal}

El modo más simple de interpolación es conectar dos puntos con una línea recta. Esta técnica, llamada interpolación lineal, la podemos representar por la siguiente formulación, la cual se calcula considerando que se tienen dos puntos $P_1= [x_1, f(x_1)]^T$ y $P_0 = [x_0, f(x_0)]^T$

\begin{equation*}
\frac{f(x)- f(x_0)}{x - x_0} = \frac{f(x_1) - f(x_0))}{x_1 - x_0} 
\end{equation*}

La cual da lugar a la siguiente ecuación

\begin{equation*}
f(x) = f(x_0) + \frac{f(x_1)-f(x_0)}{x_1 - x_0}(x - x_0)
\end{equation*}

La implementación para una función de interpolación lineal para un conjunto de puntos es:

\begin{verbatim}
function y = Interpolacion_Lineal(x0, x, y)

N =length(x);

for k=1:N-1
    if x0 >= x(k) && x0 <= x(k+1) break;
    end
end

y = y(k) + (y(k+1) - y(k))/(x(k+1)-x(k))*(x0-x(k));

end
\end{verbatim}

\subsection{Ejemplo 1} 

Dados los puntos

\begin{center}
\begin{tabular}{r|r}
$x$ & $y$ \\
\hline
0.00 & 2.10\\
1.00 & 7.70\\
2.00 & 13.60\\
3.00 & 27.20\\
4.00 & 40.90\\
5.00 & 61.10\\
\end{tabular}
\end{center}

Implementar el código correspondiente para graficar los puntos intermedios a los valores dados utilizando interpolación lineal

\begin{verbatim}
x = [0,1,2,3,4,5];
y = [2.1, 7.70, 13.60, 27.20, 40.90, 61.10];

x_nva = 0:0.01:5;

N = length(x_nva);
y_nva = zeros(N, 1);

for k=1:N
    y_nva(k) = Interpolacion_Lineal(x_nva(k), x, y);
end;

plot(x_nva, y_nva, 'r-.', x, y, 'k*');
xlabel('x');
ylabel('f(x)');
\end{verbatim}

En la figura \ref{fig:ej_int_lineal}, se muestran los puntos correspondientes calculados de acuerdo con el código implementado

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_6/figuras/ej_int_lineal.eps, width=10.0cm}}
   \caption{Un conjunto de puntos y su interpolación lineal}
   \label{fig:ej_int_lineal}
\end{figure}

\section{Interpolación cuadrática}

El error que se observa en la figura  \ref{fig:ej_int_lineal}, se debe a que suponemos que los puntos se unen por lineas. Una manera resolver este problema es suponer una función polinomial cuadrática tomando al menos tres puntos para tal propósito.  Una forma de introducir tres puntos en una función cuadrática es mediante la siguiente ecuación

\begin{equation*}
q(x) = b_0 + b_1(x - x_0) + b_2 (x - x_0)(x - x_1)
\end{equation*}

Para calcular los valores $b_0$, $b_1$ y $b_2$ necesitamos introducir información de tres puntos $[x_0,  f(x_0)]$, $[x_1, f(x_1)]$ y $[x_2, f(x_2)]$. Si sustituimos en primer $x_0$ tenemos que

\begin{equation*}
f(x_0) = b_0 + b_1(x_0 - x_0) + b_2 (x_0 - x_0)(x - x_1)
\end{equation*}

por lo tanto 
\begin{equation*}
b_0 = f(x_0)
\end{equation*}

Para calcular $b_1$ utilizamos el segundo punto así tenemos

\begin{equation*}
f(x_1) = b_0 + b_1(x_1 - x_0) + b_2 (x_1 - x_0)(x_1 - x_1)
\end{equation*}

despejando tenemos

\begin{equation*}
b_1 = \frac{ f(x_1) - f(x_0 )} {x_1 - x_0}
\end{equation*}

Para el calculo de $b_2$ sustituimos los valores de $b_0$, $b_1$ y $[x_2, f(x_2)]$ de la siguiente manera

\begin{equation*}
f(x_2) = f(x_0) + \frac{ f(x_1) - f(x_0 )} {x_1 - x_0} (x_2 - x_0) + b_2 (x_2 - x_0)(x_2 - x_1)
\end{equation*}

despejando a $b_2$

\begin{equation*}
b_2 = \frac{ f(x_2) - f(x_0) - \frac{f(x_1) - f(x_0)}{x_1 - x_0} (x_2 - x_0)}{(x_2 - x_0)(x_2 - x_1)}
\end{equation*}


\begin{equation*}
b_2 = \frac{ f(x_2) - f(x_1)  + f(x_1) - f(x_0) - \frac{f(x_1) - f(x_0)}{x_1 - x_0} (x_2 - x_0)}{(x_2 - x_0)(x_2 - x_1)}
\end{equation*}

\begin{equation*}
b_2 = \frac{ f(x_2) - f(x_1) + ( f(x_1) - f(x_0)) \left( 1 - \frac{x_2- x_0}{x_1 - x_0} \right) }{(x_2 - x_0)(x_2 - x_1)}
\end{equation*}

\begin{equation*}
b_2 = \frac{ f(x_2) - f(x_1) - ( f(x_1) - f(x_0)) \left( \frac{x_2- x_1}{x_1 - x_0} \right) }{(x_2 - x_0)(x_2 - x_1)}
\end{equation*}

finalmente obtenemos

\begin{equation*}
b_2 = \frac{\frac{f(x_2) - f(x_1)}{x_2 - x_1} - \frac{f(x_1) - f(x_0)}{x_1 - x_0}}{x_2 - x_0}
\end{equation*}

En general podemos extender el método para $N$ puntos haciendo

\[
\begin{array}{|c|c|c|c|c|c|}
  \hline
  &&&&&\\
  k  & x      &  f(x)      & b_0^{(k)}     & b_1^{(k)}                                            & b_2^{(k)}\\
  &&&&&\\
  \hline
  &&&&&\\
  0  &  x_0 &  f(x_0)  & f(x_0)  & \frac{f(x_1) - f(x_0)}{x_1 - x_0} & \frac{\frac{f(x_2) - f(x_1)}{x_2 - x_1} -  \frac{f(x_1) - f(x_0)}{x_1 - x_0}}{x_2 - x_0} \\
  &&&&&\\
  \hline
  &&&&&\\
  1  &  x_1 &  f(x_1)  & f(x_1)  & \frac{f(x_2) - f(x_1)}{x_2 - x_1} & \frac{\frac{f(x_3) - f(x_2)}{x_3 - x_2} -  \frac{f(x_2) - f(x_1)}{x_2 - x_1}}{x_3 - x_1} \\
  &&&&&\\
  \hline
  &&&&&\\
  2  &  x_2 &  f(x_2)  & f(x_2)  & \frac{f(x_3) - f(x_2)}{x_3 - x_2} & \frac{\frac{f(x_4) - f(x_3)}{x_4 - x_3} -  \frac{f(x_3) - f(x_2)}{x_3 - x_2}}{x_4 - x_2} \\
  &&&&&\\
  \hline
  &&&&&\\
  3  &  x_3 &  f(x_3)  & f(x_3)  & \frac{f(x_4) - f(x_3)}{x_4 - x_3} & \\ 
  &&&&&\\
  \hline
  &&&&&\\
  &&&&&\\
  4  &  x_4 &  f(x_4)  & f(x_4)  &   & \\ 
  &&&&&\\
  \hline
\end{array}
\]

La función quadratica correspondiente al $k$-esimo intervalo queda representada por

\begin{equation*}
q^{(k)}(x) = b_0^{(k)} + b_1^{(k)}(x-x_0) + b_2^{(k)} (x-x_1)(x-x_2)  
\end{equation*}

\subsection{Ejemplo 1}

Calcular el logaritmo de $x = 2$ dados tres puntos 

\begin{center}
\begin{tabular}{ll}
$x_0 = 1$ & $f(x_0) = 0$\\
$x_1 = 4$ & $f(x_1) = 1.386294$\\
$x_2 = 6$ & $f(x_2) = 1.791759$
\end{tabular}
\end{center}

Para resolver tenemos

\begin{equation*}
b_0 = 0
\end{equation*}

\begin{equation*}
b_1 = \frac{ f(x_1) - f(x_0 )} {x_1 - x_0} = \frac{ 1.386294 -0} {4 - 1} = 0.4620981
\end{equation*}

El coeficiente $b_2$ tenemos

\begin{equation*}
b_2 = \frac{\frac{f(x_2) - f(x_1)}{x_2 - x_1} - \frac{f(x_1) - f(x_0)}{x_1 - x_0}}{x_2 - x_0} 
\end{equation*}

\begin{equation*}
b_2 = \frac{\frac{1.791759 - 1.386294}{6 - 4} - 0.4620981}{6 - 1} = -0.0518731 
\end{equation*}

En Forma Tabular

\[
\begin{array}{|c|c|c|c|c|c|}
  \hline
  k  & x      &  f(x)      & b_0^{(k)}     & b_1^{(k)}  & b_2^{(k)}\\
  \hline
  0 & 1  & 0.000000 & 0.000000 & 0.4620980 & -0.0518731\\
  1 & 4  & 1.386294 & 1.386294 & 0.2027325 & \\
  2 & 6  & 1.791759 & 1.791759 &                   & \\
\hline
\end{array}
\]

Finalmente para calcular el valor hacemos

\begin{equation*}
q(x) = b_0 + b_1(x - x_0) + b_2 (x - x_0)(x - x_1)
\end{equation*}

\begin{equation*}
q(x) = 0 + 0.4620981(x - 1) - 0.0518731(x -1)(x - 4)
\end{equation*}

\begin{equation*}
q(2) = 0 + 0.4620981(2 - 1) - 0.0518731(2 -1)(2 - 4) = 0.5658444
\end{equation*}

En Matlab dar 

\begin{verbatim}
Interpolacion_Newton(2, [1, 4, 6] , [0, 1.386294, 1.791759])

ans =

    0.5658
\end{verbatim}

\subsection{Ejemplo 2} 

Dados los puntos

\begin{center}
\begin{tabular}{r|r}
$x$ & $y$ \\
\hline
0.00 & 2.10\\
1.00 & 7.70\\
2.00 & 13.60\\
3.00 & 27.20\\
4.00 & 40.90\\
5.00 & 61.10\\
\end{tabular}
\end{center}


En Forma Tabular los coeficientes para interpolar son:

\[
\begin{array}{|c|c|c|c|c|c|}
  \hline
  k  & x      &  f(x)      & b_0^{(k)}     & b_1^{(k)}  & b_2^{(k)}\\
  \hline
  0 &   0.00 & 2.1000   & 2.1000   &  5.6000   &  0.1500  \\
  1 &   1.00 & 7.7000   & 7.7000   &  5.9000   &  3.8500  \\
  2 &   2.00 &13.6000  & 13.6000  &  13.6000 &  0.0500  \\
  3 &   3.00 & 27.2000  & 27.2000  &  13.7000 &  3.2500 \\
  4 &   4.00  & 40.9000  & 40.9000  &  20.2000 &             \\
  5 &   5.00 & 61.1000  & 61.1000  &                &              \\
\hline
\end{array}
\]

Las funciones cuadráticas para cada intervalo quedan:

\begin{eqnarray*}
q^{(0)}(x) &=&   2.1000 +  5.6000 (x - 0) + 0.1500 (x-0)(x-1)\\
q^{(1)}(x) &=&   7.7000 +  5.9000 (x - 1) + 3.8500 (x-1)(x-2)\\
q^{(2)}(x) &=& 13.6000 + 13.6000 (x - 2) + 0.0500 (x-2)(x-3)\\
q^{(3)}(x) &=& 27.2000 + 13.7000 (x - 3) + 3.2500 (x-3)(x-4)\\
\end{eqnarray*}


Hacer una interpolación cuadrática de los datos y comparar con la interpolación lineal

\begin{verbatim}
x = [0,1,2,3,4,5];
y = [2.1, 7.70, 13.60, 27.20, 40.90, 61.10];

x_nva = 0:0.01:5;

grado = 2;

N = length(x_nva);
M = length(x);
y_nva = zeros(N, 1);

for n=1:N
    for m=1:M-1
        if x_nva(n) >= x(m) && x_nva(n) <= x(m+1) break;
        end;
    end;
    l = m+grado;
    if l > M
        l = M;
    end;
    

    y_nva(n) = Interpolacion_Newton(x_nva(n), x(m:l), y(m:l));
end;

plot(x_nva, y_nva, 'r-.', x, y, 'k*');
xlabel('x');
ylabel('f(x)');

\end{verbatim}

En la figura \ref{fig:ej_int_cuadratica}, se muestran los puntos correspondientes calculados de acuerdo con el código implementado

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_6/figuras/ej_int_cuadratica.eps, width=10.0cm}}
   \caption{Un conjunto de puntos y su interpolación cuadrática}
   \label{fig:ej_int_cuadratica}
\end{figure}

\section{Formulas de interpolación de Newton}

Las formulaciones anteriores pueden ser generalizadas para ajustar un polinomio de $n$-ésimo orden a  $n+1$ datos. El polinomio de $n$-ésimo orden es

\begin{equation}
f_n(x) = b_0 + b_1(x-x_0) + b_2(x-x0)(x-x_1) + \cdots + b_n(x-x_0)(x - x_1)\cdots(x - x_{n-1}) \label{eq:int_Newton}
\end{equation}

Como se hizo antes con las interpolaciones lineales y cuadráticas, los puntos de los datos evaluaban los coeficientes $b_0, b_1, \cdots, b_n$. Para un polinomio de $n$-ésimo orden se requiere $n+1$ puntos: $[x_0, f(x_0)]$, $[x_1, f(x_1)]$, ..., $[x_n, f(x_n)]$. Usamos estos datos y las siguientes ecuaciones para evaluar los coeficientes

\begin{eqnarray*}
b_0 &= &f(x_0)\\
b_1 &= &f[x1, x_0]\\
b_2 &= &f[x_2, x_1, x_0]\\
&\vdots &\\
b_n &= &f[x_n, x_{n-1}, \cdots, x_1, x_0]
\end{eqnarray*}

donde las evaluaciones de la función puestas entre paréntesis son diferencias divididas finitas. Por ejemplo, la primera diferencia dividida finita se representa por lo general como

\begin{equation*}
f[x_i, x_j] = \frac{f(x_i) - f(x_j)}{x_i - x_j}
\end{equation*}

La segunda diferencia dividida finita, la cual representa la diferencia de las dos primeras diferencias divididas, se expresa por lo general como

\begin{equation*}
f[x_i, x_j, x_k] = \frac{f[x_i, x_j] - f[x_j, x_k]}{x_i - x_k}
\end{equation*}

En forma similar, la $n$-ésima diferencia dividida finja es

\begin{equation*}
f[x_n, x_{n-1}, \cdots, x_1, x_0] = \frac{f[x_n, x_{n-1}, \cdots, x_1] - f[x_{n-1}, x_{n-2}, \cdots, x_1, x_0]}{x_n - x_0}
\end{equation*}

Estas diferencias pueden usarse para evaluar los coeficientes en la ecuación \ref{eq:int_Newton}, las cuales entonces se sustituirán en la siguiente ecuación, para así obtener el polinomio de interpolación

\begin{eqnarray*}
f_n(x) &=& f(x_0)+ f[x_1,x_0](x-x_0) + f[x_2, x_1, x_0] (x-x0)(x-x_1) + \\
& \cdots &  + f[x_n, x_{n-1}, \cdots, x_1, x_0] (x-x_0)(x - x_1)\cdots(x - x_{n-1}) 
\end{eqnarray*}

En la siguiente tabla se muestra el procedimiento para el calculo de las diferencias finitas divididas.

\begin{center}
\begin{tabular}{rrrrrr}
\hline
$i$ & $x_i$ & $f(x_i)$ & Primero & Segundo & Tercero\\
      &           & $b_0$  & $b_1$    &  $b_2$    &  $b_3$\\
\hline
0 & $x_0$ & $f(x_0)$ & $f[x_1, x_0]$ & $f[x_2, x_1, x_0]$ & $f[x_3, x_2, x_1, x_0]$ \\
1 & $x_1$ & $f(x_1)$ & $f[x_2, x_1]$ & $f[x_3, x_2, x_1]$ &\\
2 & $x_2$ & $f(x_2)$ & $f[x_3, x_2]$ &\\
3 & $x_3$ & $f(x_3)$ &                       &\\
\end{tabular}
\end{center}

La implementación de este algoritmo en Matlab es:

\begin{verbatim}
function yint = Interpolacion_Newton(xi, x, y) 
N = length(x);

f = zeros(N, N);

for n=1:N
    f(n, 1) = y(n);
end;

for m = 2:N
    for n = 1:N+1-m
        f(n,m)=(f(n+1,m-1)-f(n, m-1))/(x(n+m-1) - x(n));
    end;
end;


yint = f(1,1);
dx = 1;

for n=2:N
    dx = dx*(xi-x(n-1));
    yint = yint + f(1,n)*dx;
end;

end
\end{verbatim}

\subsection{Ejemplo 1}

Calcular el logaritmo de 2 utilizando la interpolación de Newton y los valores $[1, 0]$, $[4, 1.386294]$, $[6, 1.791759]$ y $[5, 1.609438]$

La solución la llevaremos a cabo de forma tabular

\begin{center}
\begin{tabular}{rrrrrr}
\hline
$i$ & $x_i$ & $f(x_i)$ & Primero & Segundo & Tercero\\
      &          &  $b_0 $ & $b_1$    & $b_2$     & $b_3$\\
\hline
    0 & 1 & 0.0000  &  0.4621  & -0.0519   &  0.0079\\
    1 & 4 & 1.3863  &  0.2027  & -0.0204       \\
    2 & 6 & 1.7918  &  0.1823  &      \\
    3 & 5 & 1.6094  &      \\
\end{tabular}
\end{center}


Para 4 puntos el polinomio de Newton de tercer orden es

\begin{equation*}
f_3(x) = 0 + 0.462098(x-1) - 0.05187311 (x-1)(x-4) + 0.00786529(x-1)(x-4)(x-6)
\end{equation*}

y el valor del logaritmo es $f_3(2) = 0.6287686$

\subsection{Ejemplo 2} 

Dados los puntos

\begin{center}
\begin{tabular}{r|r}
$x$ & $y$ \\
\hline
0.00 & 2.10\\
1.00 & 7.70\\
2.00 & 13.60\\
3.00 & 27.20\\
4.00 & 40.90\\
5.00 & 61.10\\
\end{tabular}
\end{center}

Hacer una interpolación utilizando polinomios de Newton de quinto grado de los datos y comparar con la interpolación lineal y cuadrática

\begin{verbatim}
x = [0,1,2,3,4,5];
y = [2.1, 7.70, 13.60, 27.20, 40.90, 61.10];

x_nva = 0:0.01:5;

grado = 5;

N = length(x_nva);
M = length(x);
y_nva = zeros(N, 1);

for n=1:N
    for m=1:M-1
        if x_nva(n) >= x(m) && x_nva(n) <= x(m+1) break;
        end;
    end;
    l = m+grado;
    if l > M
        l = M;
    end;
    

    y_nva(n) = Interpolacion_Newton(x_nva(n), x(m:l), y(m:l));
end;

plot(x_nva, y_nva, 'r-.', x, y, 'k*');
xlabel('x');
ylabel('f(x)');
\end{verbatim}

En la figura \ref{fig:ej_int_quinta}, se muestran los puntos correspondientes calculados de acuerdo con el código implementado

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_6/figuras/ej_int_quinta.eps, width=10.0cm}}
   \caption{Un conjunto de puntos y su interpolación con Polinomios de Newton de Quinto grado}
   \label{fig:ej_int_quinta}
\end{figure}

\section{Interpolación de Polinomios de Lagrange}

La interpolación de polinomios de Lagrange es simplemente una reformulación del polinomio de Newton que evita el cálculo por diferencias divididas. Se puede expresar de manera concisa como

\begin{equation*}
f_n(x) = \sum_{i=0}^n L_i(x) f(x_i)
\end{equation*}

donde
    
\begin{equation*}
L_i(x) = \prod_{
\begin{array}{c}
  j=0   \\
  j \neq i
\end{array}} ^n \frac{x - x_j}{x_i - x_j}
\end{equation*}

Por ejemplo la versión lineal $n=1$ es

\begin{equation*}
f_1(x) = \frac{x-x_1}{x_0 - x_1} f(x_0) + \frac{x-x_0}{x_1 - x_0} f(x_1)
\end{equation*}

y la versión de segundo orden es

\begin{equation*}
f_2(x) = \frac{(x-x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)} f(x_0) + \frac{(x-x_0)(x-x_2)}{(x_1 - x_0)(x_1 - x_2)} f(x_1) + \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} f(x_2)
\end{equation*}

El código en Matlab para este algoritmo es

\begin{verbatim}
function yint = Interpolacion_Lagrange(xint, x, y) 

suma = 0;
N = length(x);
for n=1:N
    producto = y(n);
    for m = 1:N
        if n~=m 
            producto = producto*(xint -x(m))/(x(n) - x(m));
        end;
    end;
    suma = suma + producto;
end;
yint = suma;
end
\end{verbatim}


\subsection{Ejemplo 1}

Use una interpolación del polinomio de Lagrange de primer y segundo orden para evaluar ln(2) con base en los siguientes datos $[1, 0]$, $[4, 1.386294]$ y $[6, 1.791759]$ 

La solución lineal $n=1$ es

\begin{equation*}
f_1(x) = \frac{x-4}{1 - 4} 0 + \frac{x-1}{4 -1} 1.386294 = -0.462098 + 0.462098 x
\end{equation*}

\begin{equation*}
f_1(2) = \frac{2-4}{1 - 4} 0 + \frac{2-1}{4 -1} 1.386294 = 0.4620981
\end{equation*}

y la solución de segundo $n=2$ orden es

\begin{equation*}
f_2(x) = \frac{(x-4)(x - 6)}{(1 - 4)(1 -6)} 0 + \frac{(x-1)(x-6)}{(4 -1)(4 - 6)} 1.386294 + \frac{(x-1)(x-4)}{(6-1)(6-4)} 1.791760 
\end{equation*}

La ecuación resultante de segundo grado es

\begin{equation*}
f_2(x) = -0.66959 + 0.721463 x - 0.0518731 x^2
\end{equation*}

Para calcular el valor de $ln(2)$ resolvemos
\begin{equation*}
f_2(2) = \frac{(2-4)(2 - 6)}{(1 - 4)(1 -6)} 0 + \frac{(2-1)(2-6)}{(4 -1)(4 - 6)} 1.386294 + \frac{(2-1)(2-4)}{(6-1)(6-4)} 1.791760 = 0.5658444  
\end{equation*}

Como era de esperar, ambos resultados concuerdan con los que se obtuvieron utilizando interpolación utilizando polinomios de Newton.

\subsection{Ejemplo 2}

Dados los vectores $x = [1, 2]$ y $y = [2, 4]$ calcular la interpolación de grado $n=1$.
\[
f_1(x) = L_0(x) y_0  + L_1(x) y_1
\]

\[
f_1(x) = 2 L_0(x)  + 4 L_1(x) 
\]

Los multiplicadores $L_i(x)$ se calculan con:

\[
L_0(x) = \frac{x - 2}{1 - 2} = 2-x
\]

\[
L_1(x) = \frac{x - 1}{2 - 1} = x-1
\]

Sustituyendo en la función tenemos

\[
f_1(x) = 2 (2 - x)  + 4 (x-1) = 2x  
\]

\subsection{Ejemplo 3}

Dados los puntos $x = [1, 2, 3]$  y sus correspondientes valores $y = [3, 5, 10]$ determinar:

a) El polinómio de Lagrange de grado $n = 2$

\[
f_2(x) = L_0(x) y_0 + L_1(x) y_1 + L_2(x) y_2
\]

\[
f_2(x) = 3 L_0(x) + 5 L_1(x) + 10 L_2(x) 
\]

Las funciones $L_i(x)$ son:

\[
L_0(x) \frac{(x - 2) (x - 3)}{ (1 - 2) (1 - 3)} = \frac{x^2 - 5x + 6}{2}      
\]
\[
L_1(x) = \frac{(x - 1) (x - 3)}{(2 - 1)  (2 - 3) } = -x^2 + 4x - 3
\]
\[
L_2(x) = \frac{(x - 1) (x - 2)}{(3 - 1) ( 3 - 2)}  = \frac{x^2 - 3x + 2}{2}
\]
Sustituyendo tenemos
\[
f_2(x) = 3 (x^2 - 5x + 6)/2 - 5 (x^2 - 4x +3) + 10 (x^2 - 3x + 2)/2
\]
\[
f_2(x) = 1.5x^2 -2.5x+ 4
\]

b) Corroborar que para cada valor de $x_i$ la función $f_2(x_i) = y_i$
\[
f_2(1) =  1.5 (1)^2 - 2.5 (1) + 4 = 3
\]
\[
f_2(2) = 1.5 (2)^2 - 2.5 (2) + 4 = 5
\]
\[
f_2(3) = 1.5 (3)^2 - 2.5 (3) + 4  = 10 
\]

c) Calcular el polinomio de grado 2 utilizando mínimos cuadrados para los puntos dados.

\[
M = 
\left[
\begin{array}{ccc}
   1 & 1  & 1\\
   1 & 2  & 4\\
   1 & 3  & 9
\end{array}
\right ]
\]

La solución es 

\[
a = (M^T M)^{-1} M^T y
\]

\[
a = 
\left[
\begin{array}{ccc}
  3 &  6  & 14\\
  6 & 14 & 36\\
14 & 36 & 98
\end{array}
\right ] ^{-1}
\left[
\begin{array}{c}
  18\\
  43\\
  113
\end {array}
\right ] =
\left[
\begin{array}{c}
  4\\
  -2.5\\
  1.5
\end {array}
\right ]
\]

El polinomio es :

\[
f_2(x) = 1.5 x^2 - 2.5 x + 4 
\]

En este caso la función calculada por mínimos cuadrados e interpolación de la Lagrange son exactamente iguales.

\subsection{Ejemplo 4}

Para los valores de $x = [0.5, 1.3, 2.1, 2.9]$ y $y = [1.25, 4.3, 3.9, 2.1]$ determinar:

a) La aproximación polinomial de grado $n=3$

\[
f_3(x) = 1.25 L_0(x) + 4.3 L_1(x) + 3.9 L_2(x) + 2.1 L_3(x)
\]
Las funciones $L_i(x)$ son:
\[
L_0 = \frac{(x - 1.3)(x -2.1)(x -2.9)}{(0.5 - 1.3)(0.5 - 2.1)(0.5 -2.9)} = -0.325521x^3 + 2.050781 x^2 - 4.09831 x + 2.577148  
\]

\[
L_1 = \frac{(x - 0.5)(x -2.1)(x -2.9)}{(1.3 - 0.5)(1.3 - 2.1)(1.3 -2.9)}  =  0.976563 x^3 - 5.37109 x^2 + 8.38867 x - 2.97363
\]

\[
L_2 = \frac{(x - 0.5)(x -1.3)(x -2.9)}{(2.1 - 0.5)(2.1 - 1.3)(2.1 -2.9)}  =  -0.976563x^3 + 4.58984x^2 - 5.73242 x + 1.84082
\]

\[
L_3 = \frac{(x - 0.5)(x -1.3)(x -2.1)}{(2.9 - 0.5)(2.9 - 1.3)(2.9 -2.1)}   0.3255221 x^3 - 1.26953 x^2 + 1.44206 x - 0.444336
\]

Sustituyendo tenemos

\[
f_3(x) =  
\begin{array}{l}
        1.25 (-0.325521x^3 + 2.050781x^2 - 4.09831 x + 2.577148) \\
       +4.3  (0.976563 x^3 - 5.37109 x^2 + 8.38867 x - 2.97363) \\
        + 3.9 (-0.976563x^3 + 4.58984 x^2 - 5.73242 x + 1.84082)\\
        + 2.1 (0.3255221x^3 - 1.26953 x^2 + 1.44206 x - 0.444336)
\end{array}
\]

\[
f_3(x) = 0.667321 x^3 - 5.297847 x^2 + 11.620281 x - 3.319081
\]

b) Determinar el polinomio de grado 3 de mínimos cuadrados

\[
M = 
\left [
\begin{array}{cccc}
    1.0000  &  0.5000  &  0.2500  &  0.1250\\
    1.0000  &  1.3000  &  1.6900  &  2.1970\\
    1.0000  &  2.1000  &  4.4100  &  9.2610\\
    1.0000  &  2.9000  &  8.4100  & 24.3890
\end{array}
\right ]
\]

\[ a = 
\left [
\begin{array}{cccc}
    4.0000  &  6.8000  & 14.7600 &  35.9720\\
    6.8000  & 14.7600 &  35.9720 &  93.0948\\
   14.7600 &  35.9720 &  93.0948 & 249.6967\\
   35.9720 &  93.0948 & 249.6967 & 685.4319
\end{array}
\right ]^{-1}
\left [
\begin{array}{c}
   11.5500\\
   20.4950\\
   42.4395\\
   96.9382
\end{array}
\right ] = 
\left [
\begin{array}{c}
   -3.3191 \\
   11.6203 \\
   -5.2979\\
    0.6673
\end{array}
\right ]
\]

\[
f_3(x) = 0.6673 x^3 - 5.2979 x^2 + 11.6203 - 3.3191
\]

c) Graficar la curva calculada con los puntos dados.

La Figura \ref{fig:inter_Lagrange} se muestra la grafica de la función $f_3(x)$ y los valores de $x$ y $y$. Los puntos dados se muestran como esferas de color azul y la curva calculada en rojo.

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_6/figuras/Lagrange_3.eps, width=10.0cm}}
   \caption{Aproximación de tercer grado para cuatro puntos utilizando interpolación de Lagrange}
   \label{fig:inter_Lagrange}
\end{figure}

\subsection{Ejemplo 5}

Dados los valores en la siguiente tabla calcular:
\begin{center}
\begin{tabular}{r|r}
$x$ & $f(x)$ \\
\hline
1 & 0.0000\\
2 & 0.6931\\
3 & 1.0986\\
4 & 1.3863\\
\end{tabular}
\end{center}

\begin{itemize}
\item El polinomio de Lagrange de cúbico en el intervalo $1 < x < 2$
\item El valor de interpolación en $x = 1.4$
\end{itemize}

Para calcular el Polinomio de Lagrange hacemos
\begin{eqnarray*}
f_3(x) & = & \frac{(x - 2)(x - 3)(x - 4)}{(1 - 2)(1 - 3)(1 -4)} 0.0000 +  \frac{(x - 1)(x - 3)(x - 4)}{(2 - 1)(2 - 3)(2 -4)} 0.6931 + \\
          &    & \frac{(x - 1)(x - 2)(x - 4)}{(3 - 1)(3 - 2)(3 -4)} 1.0986 + \frac{(x - 1)(x - 2)(x - 3)}{(4 - 1)(4 - 2)(4 - 3)} 1.3863 
\end{eqnarray*}

Simplificando la expresión tenemos

\begin{equation*}
f_3(x) = -1.1505 + 1.4358 x - 0.3136 x^2 + 0.0283 x^3
\end{equation*}

Para calcular el valor interpolado hacemos

\begin{eqnarray*}
f_3(1.4) & = & \frac{(1.4 - 2)(1.4 - 3)(1.4 - 4)}{(1 - 2)(1 - 3)(1 -4)} 0.0000 +  \frac{(1.4 - 1)(1.4 - 3)(1.4 - 4)}{(2 - 1)(2 - 3)(2 -4)} 0.6931 + \\
          &    & \frac{(1.4 - 1)(1.4 - 2)(1.4 - 4)}{(3 - 1)(3 - 2)(3 -4)} 1.0986 + \frac{(1.4 - 1)(1.4 - 2)(1.4 - 3)}{(4 - 1)(4 - 2)(4 - 3)} 1.3863  
\end{eqnarray*}

obviamente que la otra alternativa para evaluar es:

\begin{equation*}
f_3(1.4) = -1.1505 + 1.4358 (1.4) - 0.3136 (1.4)^2 + 0.0283 (1.4)^3
\end{equation*}

El resultado de la interpolación es:

\begin{equation*}
f_3(1.4)= 0.322619
\end{equation*}