Notebook

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.

Integración numerica¶

Introducción¶

Discutiremos de manera breve la definición de cuadratura. Posteriormente nos concentraremos en las cuadraturas Gaussianas que por su eficiencia y facilidad de sistematización son de amplio uso en ingeniería y física. Para estas cubriremos su desarrollo general y su implementación en Python. Los detalles de la cuadratura Gaussiana y su implementación se discutirán por medio de un ejemplo.

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

Identificar una cuadratura como una formula de evaluar integrales numéricamente.
Identificar la relación entre la función a integrar y el tipo de esquema requerido para su evaluación.
Evaluar integrales numéricamente usando cuadraturas Gaussianas.

Cuadraturas¶

Una cuadratura es una formula para la evaluación numerica de integrales de la forma general:

$I=\int\limits_{V(\vec{x})} f(\vec{x}) \mathrm{d}V(\vec{x}) \approx\sum_{i=1}^{N} f(\vec{x}_i)w_i\, .$

Note que esta expresión corresponde a la evaluación de la función $f(x)$ en $N$ puntos de coordenadas $x_i$ multiplicados por $N$ factores $w_i$ . Los factores se denominan pesos o factores de ponderación ya que se encargan de ponderar la contribución de cada término $f(x_i)$ a $I$ y tienen una interpretación similar al diferencial $\mathrm{d}V$ . Incluso, estos últimos son los que se encargarían de aportar las unidades pertinentes a la integral (aproximada).

Ejemplo: regla del trapecio¶

Una cuadratura con la cual estamos familiarizados es la regla del trapecio dada por:

$I=\int\limits_a^b f(x) \mathrm{d}x \approx \frac{h}{2}[f(a) + f(b)]\, ,$

en donde $h = b - a$ . En esta expresión podemos reconocer los factores de ponderación $w_1 = h/2$ , $w_2 = h/2$ y los puntos de evaluación $x_1 = a$ y $x_2 = b$ .

Por ejemplo, consideremos la siguiente integral:

$I = \int\limits_{-1}^{+1} (x^3 + 4x^2 - 10) \mathrm{d}x \approx 1.0\cdot f(-1) + 1.0\cdot f(+1) = -12\, .$

Cuadraturas Gaussianas¶

Una de las cuadraturas mas poderosas encontradas en la practica son las denominadas cuadraturas Gaussianas. En estas, los factores de ponderación $w_i$ y los puntos de evaluación $x_i$ son seleccionados de manera que se obtenga la mejor aproximación (mínimo error) de la manera más efectiva (mínimo número de puntos de evaluación). El ser formuladas usando un proceso de ajuste de $2 N$ parámetros correspondientes a los $N$ pesos y a los $N$ puntos de evaluación permiten integrar de manera exacta funciones polinomiales de orden a lo sumo $2 N - 1$ .

La principal desventaja de las cuadraturas Gaussianas es el hecho de que en estas los puntos de evaluación se encuentran especificados en términos de coordenadas en el rango fijo entre $x=-1$ y $x=+1$ lo cual obliga a que sea necesario realizar una transformación previa o cambio de variable.

Para evitar confusiones en la notación denotemos el espacio en el que se indican las coordenadas de las cuadraturas Gaussianas mediante la letra $r$ , de manera que el cambio de variables se expresa como:

$I = \int\limits_{x=a}^{x=b} f(x) \mathrm{d}x \equiv \int\limits_{r=-1}^{r=+1}F(r) \mathrm{d}r\, .$

Nótese que el cambio de variables implica:

Relacionar $x$ y $r$ lo que podemos escribir de forma general como $x = x(r)$ y $r = r(x)$ .
Expresar $f(x)$ en términos de la nueva variable de acuerdo con $F(r) = f[x(r)]$ .
Expresar $\mathrm{d}x$ en términos de $\mathrm{d}r$ .

Cuadratura de 2 puntos¶

Considere el caso de una cuadratura de 2 puntos, es decir $N =2$ . En este caso los factores de ponderación y puntos de evaluación se especifican en la siguiente tabla:

$r$	$w$
$\frac{-\sqrt3}{3}$	$1.0$
$\frac{+\sqrt3}{3}$	$1.0$

Para realizar el cambio de variables asumamos que la relación entre las variables independientes $x$ y $r$ es lineal de manera que:

$x(r) = \frac{1}{2}(a + b) + \frac{r}{2}(b - a) \equiv \frac{1}{2}(a + b) + \frac{h}{2}r\, ,$

y por lo tanto:

$\mathrm{d}x=\frac{h}{2}\mathrm{d}r\, .$

Esto que produce la siguiente equivalencia entre las integrales en los 2 espacios:

$I = \int\limits_{x=a}^{x=b} f(x) \mathrm{d}x \equiv \int\limits_{r=-1}^{r=+1} f[ x(r)]\frac{h}{2} \mathrm{d}r\, .$

Ahora, la integral formulada en el espacio de $r$ es fácilmente evaluable mediante las coordenadas y pesos de la tabla.

Consultar los factores y puntos de integración para una cuadratura Gaussiana de 4 puntos.

Solución en Python¶

En los bloques de código que se presentan a continuación se implementa la cuadratura Gaussiana de 2 puntos para calcular la integral:

$I=\int_{x = -1}^{x = +1}(x^3+4x^2-10)\operatorname dx$

Adicionar comentarios a cada uno de los bloques de código que se presentan a continuación.

In [1]:

%matplotlib notebook        
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

En el espacio encerrado entre comillas en cada una de las siguientes subrutinas indique el significado de cada uno de los parámetros y su tipo de dato.

In [2]:

def gpoints2():
    """Cuadratura de Gauss de 2 puntos"""
    xw = np.zeros([2])
    xp = np.zeros([2])
    xw[:] = 1.0
    xp[0] = -0.577350269189626
    xp[1] = 0.577350269189626
    
    return xw, xp

In [3]:

def transform(a, b, r):
    """

    """
    h = b-a
    xr = (a + b)/2.0 + h*r/2.0
        
    return xr, h

In [4]:

def myfun(x):
    """

    """
    fx = x**3 + 4*x**2 - 10
        
    return fx

Adicione comentarios al código de integración.

In [5]:

ngpts = 2
a = -1.0
b = +1.0
integral = 0.0
xw, xp = gpoints2()
for i in range(0, ngpts):
    xr, h = transform(a, b, xp[i])
    fx = myfun(xr)
    integral = integral + fx*h/2.0*xw[i]
print(integral)

-17.333333333333332

Preguntas:

Modificar el código anterior para calcular la integral con una cuadratura de 3 puntos.
Repetir el cálculo de la integral anterior si ahora los límites de integración son $a =0$ y $b=2$ .
Usando la cuadratura Gaussiana calcular la siguiente integral:

$I=\int\limits_{x=3.0}^{x=6.0} \mathrm{d}x$

¿Cómo sería la generalización de la cuadratura Gaussiana sobre un cuadrilátero?

Glosario de términos¶

Cuadratura: Formula de integración numerica compuesta por un conjunto de puntos de evaluación y factores de ponderación.

Punto de integración: Punto de evaluación de la función a integrar mediante una cuadratura numérica.

Punto de Gauss: Punto de integración en una cuadratura Gaussiana.

Factor de ponderación: Constante que pondera la contribución de la función a la integral cuando esta es evaluada en un punto de integración determinado.

Actividad para la clase¶

La figura muestra el problema de una cuña de semi-ángulo interno $\phi=\frac\pi4$ y lado $\ell = 10.0$ sometida a tracciones en las superficies inclinadas de magnitud $S = 1.0$ .

Considerando que la relaciónes deformación-desplazamiento y tensión-deformación están dadas por:

$\begin{align} \varepsilon_{xx} &= \frac{\partial u}{\partial x}\, ,\\ \varepsilon_{yy} &= \frac{\partial v}{\partial y}\, ,\\ \varepsilon_{xy} &= \frac{1}{2}\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)\, ,\\ \sigma_{xx} &= \frac E{1 + \nu}\varepsilon_{xx} + \frac{\nu E}{(1+\nu)(1-2\nu)}(\varepsilon_{xx} + \varepsilon_{yy})\, ,\\ \sigma_{yy} &= \frac E{1+\nu}\varepsilon_{yy} + \frac{\nu E}{(1+\nu)(1-2\nu)}(\varepsilon_{xx} + \varepsilon_{yy})\, ,\\ \sigma_{xy} &= \frac{E}{2(1 + \nu)} \varepsilon_{xy}\, , \end{align}$

se pide:

Calcular la energía de deformación del sistema dada por:

$I = \frac{1}{2}\int\limits_S (\sigma_{xx}\varepsilon_{xx} + \sigma_{yy}\varepsilon_{yy} + 2\sigma_{xy}\varepsilon_{xy})\mathrm{d}S\, ,$

asumiendo que los desplazamientos en los puntos izquierdo y derecho están dados por

$\vec{u}_\text{izq} = -2.0 \hat{\imath}\, ,$

$\vec{u}_\text{der} = +2.0\hat{\imath}\, ,$

mientras que los de los puntos superior e inferior corresponden a

$\vec{u}_\text{sup} = -2.0 \hat{\jmath}\, ,$

$\vec{u}_\text{inf}=+2.0\hat{\jmath}\, .$

Verifique que su resultado es correcto comparando con la solución analítica del problema.

Formato del notebook¶

La siguiente celda cambia el formato del Notebook.

In [6]:

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

Out[6]: