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}$$

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.

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.

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}$$

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.

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.

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

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$.

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

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].

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.

Paso 5: Encontrando las derivadas¶

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

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

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

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}$$

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.