Contenido bajo licencia Creative Commons BY 4.0 y código bajo licencia MIT. © Juan Gómez y Nicolás Guarín-Zapata 2020. Este material es parte del curso Modelación Computacional en el programa de Ingeniería Civil de la Universidad EAFIT.

# Interpolación de Lagrange en 1D

## Introducción

El problema de interpolación o de predecir funciones a partir de un número determinado de datos representativos de la función aparece con bastante frecuencia en la interpretación de datos experimentales. De otro lado, las técnicas de aproximación de funciones por medio de métodos de interpolación son la base para la formulación de los métodos más importantes de la mecánica computacional como lo son los elementos finitos y los elementos de frontera. En este Notebook se discutirán algunos aspectos básicos y fundamentales de la teoría de interpolación. El notebook se describe en términos del desarrollo completo de un problema ejemplo incluyendo su implementación en Python.

**Al completar este notebook usted debería estar en la capacidad de:**

* Reconocer el problema de interpolación de funciones como uno de aproximación de una función desconocida en términos de valores discretos de la misma función.

* Identificar las propiedades fundamentales de los polinomios de interpolación de Lagrange.

* Usar polinomios de interpolación de Lagrange para proponer aproximaciones a una función dado un conjunto de $N$ valores conocidos de la misma.

* Reconocer la diferencia entre variables primarias y secundarias en un esquema de interpolación.

## Interpolación de Lagrange


El problema de interpolación consiste en la determinación del valor de una función $f(x)$ en un punto arbitrario $x \in [x_1, x_n]$, dados valores conocidos de la función al interior de un dominio de solución. De acuerdo con el teorema de interpolación de Lagrange la aproximación $\hat f(x)$ a la función $f(x)$ se construye como:

\begin{equation}
\hat{f}(x) = \sum_{I=1}^n L_I(x) f_I
\end{equation}

donde $L_I$ es el $I-$ésimo polinomio de orden $n-1$ y $f_1, f_2, \cdots, f_n$ son los $n$ valores conocidos de la función. El $I-$ésimo polinomio de Lagrange se calcula siguiendo la siguiente expresión:

\begin{equation}
{L_I}(x) = \prod_{J=1, J \ne I}^{n} \frac{(x - x_J)}{(x_I - x_J)}\, .
\end{equation}

### Primera derivada

En varios problemas de ingeniería es necesario usar técnicas de interpolación para encontrar valores de las derivadas de la variable primaria o principal. Por ejemplo, en problemas de mecánica de sólidos y teoría de elasticidad, en los cuales la variable primaria es el campo de desplazamientos, uno puede estar interesado en determinar las deformaciones unitarias las cuales son derivadas espaciales de los desplazamientos. Considerando que solo se dispone de valores discretos de los desplazamientos se tiene que las derivadas que se requieren se pueden encontrar usando estos valores discretos. Lo anterior equivale a derivar $\hat{f}(x)$ directamente como sigue:


\begin{equation}
\frac{\mathrm{d}\hat{f}}{\mathrm{d}x}=\frac{\mathrm{d}L_I(x)}{\mathrm{d}x}f_I
\end{equation}

Haciendo,

$$B_I(x) = \frac{\mathrm{d}L_I(x)}{\mathrm{d}x}\, ,$$

podemos escribir el esquema de interpolación como:

\begin{equation}
\frac{\mathrm{d}\hat{f}}{\mathrm{d}x} = \sum_{I=1}^n B_I(x) f_I\, .
\end{equation}

## Ejemplo





Formule un esquema de interpolación para encontrar un valore de la función

\begin{equation}
f(x) = x^3 + 4x^2 - 10
\end{equation}

en un punto arbitrario $x$ en el intervalo $[ -1, 1]$ asumiendo que se conoce el valor exacto de la función en los puntos $x=-1.0$, $x=+1.0$ y $x=0.0.$

En este ejemplo se conoce la función y aparentemente no tiene mucho sentido buscar una aproximación de la misma resolviendo un problema de interpolación. Sin embargo es conveniente seleccionar una función conocida para poder asimilar el problema numérico y sus limitaciones. En este contexto asumiremos que en una aplicación real conocemos los valores de la función en un conjunto de puntos $x=-1.0$, $x=+1.0$ and $x=0.0$ los cuales se denominan puntos nodales o simplemente _nodos_.

El proceso de interpolación involucra 2 pasos fundamentales: 

1. Determinar los polinomios de interpolación $L_I$ usando la productoria.

2. Usar la combinación lineal para construir el polinomio interpolante o la aproximación a la función $\hat f(x)$.

Veamos estos pasos entonces.

1. Considerando que tenemos 3 puntos _nodales_ necesitamos generar 3 polinomios de interpolación de segundo orden, cada uno de ellos asociado a cada punto nodal. Rotulemos los _nodos_ como $x_0 = -1.0$, $x_1 = +1.0$ y $x_2 = 0.0$. De acuerdo con esta denominación $L_0(x)$, $L_1(x)$ y $L_2(x)$ serán los polinomios de interpolación de segundo orden asociados a los puntos nodales $x_0 = -1.0$, $x_1 = +1.0$ y $x_2 = 0.0$. Usando la fórmula de la productoria tenemos:

\begin{align}
&L_0(x) = \frac{(x - x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)} \equiv \frac{1}{2}(x - 1.0)x\\
&L_1(x) = \frac{(x-x^0)(x-x^2)}{(x^1-x^0)(x^1-x^2)}\equiv\frac12(x+1.0)x\\
&L_2(x) = \frac{(x-x^0)(x-x^1)}{(x^2-x^0)(x^2-x^1)}\equiv-(x+1.0)(x-1.0)\, .
\end{align}

2. Para llegar a la aproximación final de la función realizamos la combinación lineal:

\begin{equation}
\hat{f}(x) = L_0(x)f_0 + L_1(x)f_1 + L_2(x)f_2
\end{equation}

<div class="alert alert-warning">

**Preguntas**
    
- Verificar que los polinomios de interpolación $L_0(x)$, $L_1(x)$ y $L_2(x)$ satisfacen la propiedad

$$L_I(x_J)= \delta_{IJ} \equiv \begin{cases}
1\quad \text{si } I=J\\
0\quad \text{si } I\neq J
\end{cases}\, .$$

- Si la condición $L_I(x_J)=\delta_{IJ}$ no se satisface por una de las funciones de interpolación encontradas, ¿cuál sería el efecto resultante en la aproximación de la función?
</div>

### Algunas notas de interés


* En el método de los elementos finitos es común denominar a los polinomios de interpolación como _funciones de forma_.

* El dominio de calculo es el intervalo de tamaño 2 comprendido entre $x=-1.0$ y $x=+1.0$. Como se discutirá mas adelante, por facilidades en la programación en la implementación de métodos de elementos finitos los dominios de tamaño diferente son transformados al dominio de tamaño 2.

## Solución en Python



En los siguientes bloques de código mostramos el paso a paso en la construcción del polinomio de interpolación final $f(x)$ usando Python. Para seguir el Notebook se recomienda que de manera simultánea se implementen los bloques de código en un script independiente o que añada comentarios a las instrucciones más relevantes.

### Paso 1: Importación de módulos


En la escritura de scripts en Python es necesario importar **módulos** o bibliotecas (algunas personas usan la palabra _librería_ como mala traducción de _library_ del inglés) que contienen funciones de Python predefinidas. En este caso importaremos los **módulos**:

* `numpy` el cual es una biblioteca de funciones para realizar operaciones con matrices similar a Matlab.
* `matplotlib` el cual es una biblioteca de graficación.
* `scipy` el cual es una biblioteca fundamental para computación científica.
* `sympy` el cual es una biblioteca para realizar matemáticas simbólicas.

Python importa los módulos usando la palabra reservada `import` seguida del nombre del módulo y un nombre corto para ser usado como prefijo en referencias posteriores a las funciones contenidas en ese módulo.


```python
%matplotlib notebook
```


```python
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
import sympy as sym
```


```python
sym.init_printing()
```

### Paso 2: Definición de una función para determinar los polinomios de interpolación de Lagrange



En un programa de computador una **función** (o también llamada **subrutina**) es un bloque de código que realiza una tarea específica múltiples veces dentro de un programa o probablemente en diferentes programas. En Python estas funciones se definen por medio de la palabra clave `def` seguida del nombre dado a la función.

Conjuntamente con el nombre, encerrado entre paréntesis, la definición de la función debe incluir también una lista de parámetros (o argumentos) a ser usados por la función cuando esta realice tareas especificas.

En el ejemplo definiremos una función de Python para generar el polinomio de Lagrange usando la productoria definida previamente. Le daremos a esta función el nombre `lagrange_poly`. Sus parámetros de entrada serán la variable independiente $x$ a ser usada en el polinomio; el orden del polinomio definido por `order`; y el punto `i` asociado al polinomio. La función será usada posteriormente desde el programa principal.


```python
def lagrange_poly(x, order, i, xi=None): 
    if xi == None:
        xi = sym.symbols('x:%d'%(order+1))
    index = list(range(order+1))
    index.pop(i)
    return sym.prod([(x - xi[j])/(xi[i] - xi[j]) for j in index])
```

Alternativamente a la definición de una función usando la palabra clave `def` Python también permite la definición de funciones en el sentido del Calculo, es decir definición de relaciones que permiten  transformar escalares o vectores en otros escalares o vectores. Esto es posible a través de la opción `lambda`. En este ejemplo usaremos la opción `lambda` para crear la función:

\begin{equation}
f(x)=x^3+4x^2-10.
\end{equation}


<div class="alert alert-warning">

**Preguntas**

Use una terminal o un script independiente para probar el uso de la opción `lambda` en la definición de una función. Intente con diferentes funciones.
</div>


```python
fun = lambda x: x**3 + 4.0*x**2 - 10.0
```

Una vez creada la función $f(x)$ usando la opción `lambda` pasamos a definir un conjunto de puntos de evaluación. El numero de puntos se define por medio de la variable `npts` y usamos la función `linspace` del módulo `numpy` para crear una distribución equidistante de puntos entre $x = -1.0$ y $x = +1.0.$ Simultáneamente, el arreglo nulo `yy` almacenará los valores de la función en los puntos almacenados en `xx`.

Note que Python inicia desde la posición 0 el conteo de elementos en arreglos y otras estructuras de datos, de manera que empezamos a contar desde cero para mantener la consistencia con la implementación.


```python
npts = 200
x_pts = np.linspace(-1, 1, npts)
```

Con la función ahora disponible podemos calcular los valores (conocidos) listos para ser usados en el esquema de interpolación. Estos valores se almacenarán en el arreglo `fd()`. Para calcular cada valor de la función usamos el comando (ya disponible) `fx` correspondiente al nombre que usamos en la declaración de la función usando la opción `lambda`.


<div class="alert alert-warning">

**Preguntas**

Intente usar nombres diferentes en la definición de la función usando la opción `lambda` y ensaye si su código aún funciona.

</div>


```python
pts = np.array([-1, 1, 0])
fd = fun(pts)
```

Ahora evalúe la función en los `npts` puntos y grafique la misma:


```python
plt.figure()
y_pts = fx(x_pts)
plt.plot(x_pts , y_pts)
plt.plot(pts, fd, 'ko')
```


    <IPython.core.display.Javascript object>









    [<matplotlib.lines.Line2D at 0x7f982bdebbe0>]



### Paso 3: Encontrando los polinomios de interpolación de Lagrange



Calculemos ahora los polinomios de Lagrange invocando la función `lagrange_poly()`. Crearemos una lista vacía que llamaremos `pol` y cada que determinemos un nuevo polinomio lo adicionaremos a la lista usando para ello el método `.append()`. En ese momento los polinomios almacenados en la lista `pol=[]` existirán en notación simbólica (como si estuviéramos resolviendo el problema a mano) y estarán listos para ser evaluados en valores específicos de $x$.


```python
x = sym.symbols('x')                                       
pol = []                                                
pol.append(sym.simplify(lagrange_poly(x, 2, 0, [-1,1,0])))  
pol.append(sym.simplify(lagrange_poly(x, 2, 1, [-1,1,0])))
pol.append(sym.simplify(lagrange_poly(x, 2, 2, [-1,1,0])))
pol
```

<div class="alert alert-warning">

**Preguntas**

- Grafique los 3 polinomios de Lagrange de orden 2 y verifique que satisfacen la propiedad $L_I(x_J)=\delta_{IJ}$.

- Adicione más puntos al dominio de interpolación y grafique los polinomios resultantes. ¿Cuál es el efecto en los polinomios?
</div>

El método `subs` del módulo `sympy` substituye el valor de la variable independiente $x$ por el valor almacenado en $xx[i]$.


```python
plt.figure()
                      
for k in range(3):
    for i in range(npts):
        yy[i] = pol[k].subs([(x, xx[i])])   
    plt.plot(xx, yy)
```


    <IPython.core.display.Javascript object>






### Paso 4: Encontrando el polinomio de interpolación para aproximar la función  $f(x)$



Construyamos ahora el polinomio de aproximación completo $\hat f(x)$ de acuerdo con la superposición lineal

$$\hat{f}(x) = L_0(x)f(x_0) + L_1(x) f(x_1) +  L_2(x) f(x_2)\, ,$$


y utilizando la lista (ya disponible) `pol` que almacena los 3 polinomios generados. Solo para efectos de ilustrar primero lo mostramos con `display`. La versión evaluada de $\hat{f}(x)$ será almacenada en el arreglo `yy[i]`.


```python
display(pol[0]*fd[0] + pol[1]*fd[1] + pol[2]*fd[2])
```


```python
plt.figure()                                                                      
for i in range(npts):
    yy[i] = fd[0]*pol[0].subs([(x, xx[i])]) + fd[1]*pol[1].subs([(x, xx[i])]) \
            + fd[2]*pol[2].subs([(x, xx[i])])

zz = fx(xx)
plt.plot([-1, 1, 0], fd, 'ko')
plt.plot(xx, yy , 'r--')
plt.plot(xx, zz)
```


    <IPython.core.display.Javascript object>









    [<matplotlib.lines.Line2D at 0x7f982bea99e8>]



Note que debido a la diferencia en orden entre la aproximación $\hat{f}(x)$ y la función $f(x)$  el polinomio de interpolación no coincide completamente con la función aunque si reproduce los valores correctos en los puntos nodales.

<div class="alert alert-warning">

**Pregunta**

¿Cómo podemos mejorar la aproximación $\hat{f}(x)$  a la función $f(x)$?
</div>


### Paso 5: Encontrando las derivadas


Usemos la opción `lambda` una vez más para definir una nueva función `fdx`correspondiente a la primera derivada:

$$f'(x) = 3x^2 + 8x\, .$$

Las derivadas en los puntos nodales se almacenarán en el arreglo `fc`


```python
fdx = lambda x: 3*x**2+8.0*x
fc = np.array([fdx(-1.0), fdx(1.0) ,fdx(0.0)])
```

La interpolación de las derivadas se calculan de acuerdo con

\begin{equation}
\hat{f}'(x) = \frac{\mathrm{d}L_0(x)}{\mathrm{d}x} f_0 + \frac{\mathrm{d}L_1(x)}{\mathrm{d}x} f_1 + \frac{\mathrm{d}L_2(x)}{\mathrm{d}x} f_2\, .
\end{equation}


```python
dpol = []
dpol.append(sym.diff(pol[0],x))     
dpol.append(sym.diff(pol[1],x))
dpol.append(sym.diff(pol[2],x))
display(dpol)
display(dpol[0]*fd[0] + dpol[1]*fd[1] + dpol[2]*fd[2])

plt.figure()
yy= fdx(xx)
plt.plot(xx, yy ,'r--')
plt.plot([-1, 1, 0], fc, 'ko')

for i in range(npts):
    yy[i] = fd[0]*dpol[0].subs([(x, xx[i])]) + fd[1]*dpol[1].subs([(x, xx[i])])\
            + fd[2]*dpol[2].subs([(x, xx[i])]) 
plt.plot(xx, yy)
```

## Glosario de términos



**Fórmula de productoria:** Expresión usada para calcular los polinomios de interpolación de Lagrange.

**Función de forma:** Nombre dado a un polinomio de interpolación en el contexto del método de los elementos finitos.

**Punto nodal:** (También nodo). Nombre dado a los puntos específicos donde se conoce el valor de una función y usado en la construcción del esquema de interpolación.

**Subrutina:** Bloque de código independiente que realiza una tarea de cómputo determinada dentro de un programa principal.

**Matriz de interpolación:** Arreglo unidimensional o bidimensional que almacena las funciones de forma en un esquema de interpolación determinado.

## Actividad para la Clase

El propósito de esta actividad es familiarizar al estudiante con el uso del método de interpolación de Lagrange en un contexto particular de la ingeniería, como el método de los elementos finitos. El método de elementos finitos utiliza dominios de interpolación de tamaño constante (llamados **elementos**) permitiendo el uso de funciones de interpolación fijos y al mismo tiempo favoreciendo la automatización en programas de computador. En una aplicación típica de elementos finitos los valores nodales de una función (por ejemplo el desplazamiento) son determinados tras resolver un sistema de ecuaciones algebraicas. Estos valores nodales son posteriormente usados para encontrar el desplazamiento a lo largo de todo el dominio del problema, y también para realizar cálculos posteriores y obtener información adicional del problema físico.

**Siga los pasos que se indican a continuación para implementar, en un notebook independiente, un esquema de interpolación típico del método de elementos finitos. El esquema de interpolación resultante será el correspondiente a un solo elemento**

* Asumiendo que un dominio de interpolación constante se define como $x \in [-1.0, +1.0]$ con puntos nodales correspondientes a $x_0 = -1.0$, $x_1 = +1.0$, $x_2 = 0.0$, $x_3 = -0.25$ y $x_4 = +0.25$ use la función `lagrange_poly()` para generar las funciones de interpolación asociadas a estos 5 puntos nodales. Presente los polinomios resultantes en una gráfica.

* Utilice las funciones de interpolación encontradas en el paso anterior, para definir una matriz $[N(x)]$, que se denominará en adelante **Matriz de interpolación**, de manera que la operación de interpolar se pueda realizar usando operaciones matriciales como:

  $$\hat{f}(x) = [N(x)] \{F\}\, ,$$

  y en la cual $F$ es un vector que almacena los valores nodales de la función.

* En su Notebook imprima la versión simbólica de la matriz de interpolación.

* En su Notebook grafique una función $f(x)$ (representando el desplazamiento de una barra) y su versión aproximada $\hat{f}(x)$  usando el esquema matricial desarrollado en el numeral anterior. El código también debe graficar la primera derivada de la función.

* Para realizar la interpolación en el Notebook programe una función o subrutina de Python que reciba como parámetros de entrada las coordenadas $x$ donde se desea calcular una aproximación de la función y el vector de valores nodales de la función de desplazamientos y entregue después de su ejecución el valor interpolado de la función.

## Formato del notebook

La siguiente celda cambia el formato del Notebook.


```python
from IPython.core.display import HTML
def css_styling():
    styles = open('./nb_style.css', 'r').read()
    return HTML(styles)
css_styling()
```





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

<style>

/*
Template for Notebooks for Modelación computacional.

Based on Lorena Barba template available at:

    https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css
*/

/* Fonts */
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

/* Text */
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}	
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}	
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}
</style>








475e8a6f37d806e631e19bd108d1c1aa0abe34d1
ipynb
Jupyter Notebook
notebooks/02a_interpolacion.ipynb
AppliedMechanics-EAFIT/Mod_Temporal
6a0506d906ed42b143b773777e8dc0da5af763eb
MIT
2019-02-20T18:14:01.000Z
2020-07-19T22:44:44.000Z
notebooks/02a_interpolacion.ipynb
AppliedMechanics-EAFIT/Mod_Temporal
6a0506d906ed42b143b773777e8dc0da5af763eb
MIT
2020-04-15T00:22:58.000Z
2020-07-04T17:03:54.000Z
notebooks/02a_interpolacion.ipynb
AppliedMechanics-EAFIT/Mod_Temporal
6a0506d906ed42b143b773777e8dc0da5af763eb
MIT
2020-05-14T18:17:09.000Z
2020-10-27T06:37:05.000Z
Qwen/Qwen-72B
1. YES
2. YES