{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Integración numerica"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introducción"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"**Al completar este notebook usted debería estar en la capacidad de:**\n",
"\n",
"* Identificar una cuadratura como una formula de evaluar integrales numéricamente.\n",
"\n",
"* Identificar la relación entre la función a integrar y el tipo de esquema requerido para su evaluación.\n",
"\n",
"* Evaluar integrales numéricamente usando cuadraturas Gaussianas."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cuadraturas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Una cuadratura es una formula para la evaluación numerica de integrales de la forma general:\n",
"\n",
"\n",
"$$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\\, .$$\n",
"\n",
"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)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Ejemplo: regla del trapecio\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Una cuadratura con la cual estamos familiarizados es la regla del trapecio dada por:\n",
"\n",
"$$I=\\int\\limits_a^b f(x) \\mathrm{d}x \\approx \\frac{h}{2}[f(a) + f(b)]\\, ,$$\n",
"\n",
"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$.\n",
"\n",
"Por ejemplo, consideremos la siguiente integral:\n",
"\n",
"$$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\\, .$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Cuadraturas Gaussianas\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Una de las cuadraturas mas poderosas encontradas en la practica son las denominadas cuadraturas [Gaussianas](https://en.wikipedia.org/wiki/Gaussian_quadrature). 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$.\n",
"\n",
"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.\n",
"\n",
"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:\n",
"\n",
"$$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",
"\n",
"Nótese que el cambio de variables implica:\n",
"\n",
"* Relacionar $x$ y $r$ lo que podemos escribir de forma general como $x = x(r)$ y $r = r(x)$.\n",
"\n",
"* Expresar $f(x)$ en términos de la nueva variable de acuerdo con $F(r) = f[x(r)]$.\n",
"\n",
"* Expresar $\\mathrm{d}x$ en términos de $\\mathrm{d}r$.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cuadratura de 2 puntos"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:\n",
"\n",
"\n",
"| $r$ | $w$ |\n",
"|---------------------|-------|\n",
"| $\\frac{-\\sqrt3}{3}$ | $1.0$ |\n",
"| $\\frac{+\\sqrt3}{3}$ | $1.0$ |\n",
"\n",
"\n",
"Para realizar el cambio de variables asumamos que la relación entre las variables independientes $x$ y $r$ es lineal de manera que:\n",
"\n",
"$$x(r) = \\frac{1}{2}(a + b) + \\frac{r}{2}(b - a) \\equiv \\frac{1}{2}(a + b) + \\frac{h}{2}r\\, ,$$\n",
"\n",
"y por lo tanto:\n",
"\n",
"$$\\mathrm{d}x=\\frac{h}{2}\\mathrm{d}r\\, .$$\n",
"\n",
"Esto que produce la siguiente equivalencia entre las integrales en los 2 espacios:\n",
"\n",
"$$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\\, .$$\n",
"\n",
"Ahora, la integral formulada en el espacio de $r$ es fácilmente evaluable mediante las coordenadas y pesos de la tabla."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div class=\"alert alert-warning\">\n",
"Consultar los factores y puntos de integración para una cuadratura Gaussiana de 4 puntos.\n",
"</div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solución en Python\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"En los bloques de código que se presentan a continuación se implementa la cuadratura Gaussiana de 2 puntos para calcular la integral:\n",
"\n",
"$$\n",
"I=\\int_{x = -1}^{x = +1}(x^3+4x^2-10)\\operatorname dx\n",
"$$\n",
"\n",
"<div class=\"alert alert-warning\">\n",
"Adicionar comentarios a cada uno de los bloques de código que se presentan a continuación.\n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib notebook \n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import sympy as sym"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div class=\"alert alert-warning\">\n",
"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.\n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def gpoints2():\n",
" \"\"\"Cuadratura de Gauss de 2 puntos\"\"\"\n",
" xw = np.zeros([2])\n",
" xp = np.zeros([2])\n",
" xw[:] = 1.0\n",
" xp[0] = -0.577350269189626\n",
" xp[1] = 0.577350269189626\n",
" \n",
" return xw, xp"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def transform(a, b, r):\n",
" \"\"\"\n",
"\n",
" \"\"\"\n",
" h = b-a\n",
" xr = (a + b)/2.0 + h*r/2.0\n",
" \n",
" return xr, h"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def myfun(x):\n",
" \"\"\"\n",
"\n",
" \"\"\"\n",
" fx = x**3 + 4*x**2 - 10\n",
" \n",
" return fx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div class=\"alert alert-warning\">\n",
"Adicione comentarios al código de integración.\n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-17.333333333333332\n"
]
}
],
"source": [
"ngpts = 2\n",
"a = -1.0\n",
"b = +1.0\n",
"integral = 0.0\n",
"xw, xp = gpoints2()\n",
"for i in range(0, ngpts):\n",
" xr, h = transform(a, b, xp[i])\n",
" fx = myfun(xr)\n",
" integral = integral + fx*h/2.0*xw[i]\n",
"print(integral)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div class=\"alert alert-warning\">\n",
"\n",
"**Preguntas:**\n",
"\n",
"1. Modificar el código anterior para calcular la integral con una cuadratura de 3 puntos.\n",
"\n",
"2. Repetir el cálculo de la integral anterior si ahora los límites de integración son $a =0$ y $b=2$.\n",
"\n",
"3. Usando la cuadratura Gaussiana calcular la siguiente integral:\n",
"\n",
"$$I=\\int\\limits_{x=3.0}^{x=6.0} \\mathrm{d}x$$\n",
"\n",
"4. ¿Cómo sería la generalización de la cuadratura Gaussiana sobre un cuadrilátero?\n",
"\n",
"</div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Glosario de términos"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Cuadratura:** Formula de integración numerica compuesta por un conjunto de puntos de evaluación y factores de ponderación.\n",
"\n",
"**Punto de integración:** Punto de evaluación de la función a integrar mediante una cuadratura numérica.\n",
"\n",
"**Punto de Gauss:** Punto de integración en una cuadratura Gaussiana.\n",
"\n",
"**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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Actividad para la clase"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"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$.\n",
"\n",
"\n",
"<center><img src=\"img/wedge.png\"\n",
" alt=\"Esquema de la cuña.\"\n",
" style=\"width:300px\">\n",
"</center>\n",
"\n",
"\n",
"Considerando que la relaciónes deformación-desplazamiento y tensión-deformación están dadas por:\n",
"\n",
"\\begin{align}\n",
"\\varepsilon_{xx} &= \\frac{\\partial u}{\\partial x}\\, ,\\\\\n",
"\\varepsilon_{yy} &= \\frac{\\partial v}{\\partial y}\\, ,\\\\\n",
"\\varepsilon_{xy} &= \\frac{1}{2}\\left(\\frac{\\partial u}{\\partial y}\n",
" + \\frac{\\partial v}{\\partial x}\\right)\\, ,\\\\\n",
"\\sigma_{xx} &= \\frac E{1 + \\nu}\\varepsilon_{xx} + \\frac{\\nu E}{(1+\\nu)(1-2\\nu)}(\\varepsilon_{xx} + \\varepsilon_{yy})\\, ,\\\\\n",
"\\sigma_{yy} &= \\frac E{1+\\nu}\\varepsilon_{yy} + \\frac{\\nu E}{(1+\\nu)(1-2\\nu)}(\\varepsilon_{xx} + \\varepsilon_{yy})\\, ,\\\\\n",
"\\sigma_{xy} &= \\frac{E}{2(1 + \\nu)} \\varepsilon_{xy}\\, ,\n",
"\\end{align}\n",
"\n",
"se pide:\n",
"\n",
"1. Calcular la energía de deformación del sistema dada por:\n",
"\n",
"$$I = \\frac{1}{2}\\int\\limits_S (\\sigma_{xx}\\varepsilon_{xx} + \\sigma_{yy}\\varepsilon_{yy}\n",
" + 2\\sigma_{xy}\\varepsilon_{xy})\\mathrm{d}S\\, ,$$\n",
"\n",
"asumiendo que los desplazamientos en los puntos izquierdo y derecho están dados por \n",
"\n",
"$$\\vec{u}_\\text{izq} = -2.0 \\hat{\\imath}\\, ,$$\n",
"\n",
"y\n",
"\n",
"$$\\vec{u}_\\text{der} = +2.0\\hat{\\imath}\\, ,$$\n",
"\n",
"mientras que los de los puntos superior e inferior corresponden a\n",
"\n",
"$$\\vec{u}_\\text{sup} = -2.0 \\hat{\\jmath}\\, ,$$\n",
"\n",
"y\n",
"\n",
"$$\\vec{u}_\\text{inf}=+2.0\\hat{\\jmath}\\, .$$\n",
"\n",
"2. Verifique que su resultado es correcto comparando con la solución analítica del problema."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Formato del notebook"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La siguiente celda cambia el formato del Notebook."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
"<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'>\n",
"<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
"\n",
"<style>\n",
"\n",
"/*\n",
"Template for Notebooks for Modelación computacional.\n",
"\n",
"Based on Lorena Barba template available at:\n",
"\n",
" https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css\n",
"*/\n",
"\n",
"/* Fonts */\n",
"@font-face {\n",
"font-family: \"Computer Modern\";\n",
"src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
"}\n",
"\n",
"/* Text */\n",
"div.cell{\n",
"width:800px;\n",
"margin-left:16% !important;\n",
"margin-right:auto;\n",
"}\n",
"h1 {\n",
"font-family: 'Alegreya Sans', sans-serif;\n",
"}\n",
"h2 {\n",
"font-family: 'Fenix', serif;\n",
"}\n",
"h3{\n",
"font-family: 'Fenix', serif;\n",
"margin-top:12px;\n",
"margin-bottom: 3px;\n",
"}\n",
"h4{\n",
"font-family: 'Fenix', serif;\n",
"}\n",
"h5 {\n",
"font-family: 'Alegreya Sans', sans-serif;\n",
"}\t\n",
"div.text_cell_render{\n",
"font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
"line-height: 135%;\n",
"font-size: 120%;\n",
"width:600px;\n",
"margin-left:auto;\n",
"margin-right:auto;\n",
"}\n",
".CodeMirror{\n",
"font-family: \"Source Code Pro\";\n",
"font-size: 90%;\n",
"}\n",
"/* .prompt{\n",
"display: None;\n",
"}*/\n",
".text_cell_render h1 {\n",
"font-weight: 200;\n",
"font-size: 50pt;\n",
"line-height: 100%;\n",
"color:#CD2305;\n",
"margin-bottom: 0.5em;\n",
"margin-top: 0.5em;\n",
"display: block;\n",
"}\t\n",
".text_cell_render h5 {\n",
"font-weight: 300;\n",
"font-size: 16pt;\n",
"color: #CD2305;\n",
"font-style: italic;\n",
"margin-bottom: .5em;\n",
"margin-top: 0.5em;\n",
"display: block;\n",
"}\n",
".warning{\n",
"color: rgb( 240, 20, 20 )\n",
"}\n",
"</style>\n",
"\n",
"<script>\n",
"/* Equations */\n",
"\n",
"MathJax.Hub.Config({\n",
"TeX: {\n",
"extensions: [\"AMSmath.js\"]\n",
"},\n",
"tex2jax: {\n",
"inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
"displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
"},\n",
"displayAlign: 'center', // Change this to 'center' to center equations.\n",
"\"HTML-CSS\": {\n",
"styles: {'.MathJax_Display': {\"margin\": 4}}\n",
"}\n",
"});\n",
"</script>\n",
"\n",
"\n"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = open('./nb_style.css', 'r').read()\n",
" return HTML(styles)\n",
"css_styling()"
]
}
],
"metadata": {
"celltoolbar": "Raw Cell Format",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.10"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}