\n",
"\n",
"**Preguntas**\n",
" \n",
"- Modificar los ejes de la gráfica y adicionar una cuadricula para facilitar la identificación aproximada de las raíces.\n",
"\n",
"- ¿Cuáles son los valores aproximados de las raíces que se obtienen de la gráfica?\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Para hacer la búsqueda inicial usaremos el algoritmo de _bracketing_ el cual recibe como parámetros de entrada la función, los extremos del intervalo de búsqueda y el número de intervalos a considerar. En el siguiente bloque de código definimos la función **bracketing()** de manera que esta pueda ser utilizada de forma general."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Comentar apropiadamente el siguiente código\n",
"def bracketing(fun, a, b, N):\n",
" \"\"\"Estimate the root of a real function using bracketing\n",
"\n",
" Parameters\n",
" ----------\n",
" fun : function\n",
" Function.\n",
" a : float\n",
" Initial point of the interval.\n",
" b : float\n",
" Final point of the interval.\n",
" N : int\n",
" Number of subdivisions for the interval\n",
"\n",
" Returns\n",
" -------\n",
" xR : list\n",
" List with intervals where a change of sign was found.\n",
" msg : string\n",
" Message regarding the success of the method.\n",
" \"\"\"\n",
" msg = \"Maximum number of iterations reached.\"\n",
" dx = (b - a)/(N - 1)\n",
" iroot = 0\n",
" x2 = a\n",
" xR = []\n",
" for i in range(0, N):\n",
" x1 = x2\n",
" x2 = x1 + dx\n",
" if (fun(x1) * fun(x2)) < 0:\n",
" msg = \"A change of sign was found.\"\n",
" iroot = iroot + 1\n",
" xR.append([x1, x2])\n",
" return xR, msg"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"En el uso de la rutina de búsqueda partiremos el dominio del problema en $N = 100$ sub-intervalos de manera que tendremos:\n",
"\n",
"$$\\triangle\\sigma=\\frac{b-a}{100}\\equiv\\frac{1000}{100}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"**Preguntas**\n",
" \n",
"- ¿Qué tan cercanas son las raíces encontradas por la rutina con las identificadas en la gráfica?\n",
"\n",
"- ¿Qué pasa si se cambia el valor del numero de sub-intervalos en la operación de bracketing?\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[[-277.7777777777778, -166.6666666666667],\n",
" [-55.5555555555556, 55.555555555555515],\n",
" [388.8888888888888, 499.9999999999999]]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xR, msg = bracketing(fun, -500, 500, 10)\n",
"xR"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Paso 2a: Afinación de la raíz por el método de bisección"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"El método de bisección consiste en la partición por mitades e identificación de la mitad en la que se presenta el cambio de signo, como se muestra en la figura. El método cada vez toma un sub-intervalo de menor tamaño variando la posición del punto $x_3$.\n",
"\n",
"\n",
"
\n",
" \n",
"
\n",
"\n",
"A continuación se presenta el método de bisección que toma como parámetros de entrada la función; los valores extremos del intervalo de búsqueda $a$ y $b$; y 2 valores de tolerancia usados para detectar la presencia de una raíz ya sea por el tamaño del intervalo de búsqueda o por la cercanía con cero del valor de la función."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# Comentar apropiadamente el siguiente código\n",
"def bisection(fun, a, b, niter=50, ftol=1e-12, verbose=False):\n",
" \"\"\"Use bisection method to estimate the root of a real function\n",
"\n",
" Parameters\n",
" ----------\n",
" fun : function\n",
" Function.\n",
" grad : function\n",
" Derivative of the function.\n",
" a : float\n",
" Initial point of the interval.\n",
" b : float\n",
" Final point of the interval.\n",
" niter : int, optional\n",
" Maximum number of iterations.\n",
" ftol : float, optional\n",
" Tolerance for the root.\n",
" verbose : bool, optional\n",
" If True, prints each iteration.\n",
"\n",
" Returns\n",
" -------\n",
" x : float\n",
" Approximated root.\n",
" msg : string\n",
" Message regarding the success of the method.\n",
" \"\"\"\n",
" if fun(a) * fun(b) > 0:\n",
" c = None\n",
" msg = \"The function should have a sign change in the interval.\"\n",
" else:\n",
" for cont in range(niter):\n",
" c = 0.5*(a + b)\n",
" if verbose:\n",
" print(\"n: {}, x: {}\".format(cont, c))\n",
" if abs(fun(c)) < ftol:\n",
" msg = \"Root found with desired accuracy.\"\n",
" break\n",
" elif fun(a) * fun(c) < 0:\n",
" b = c\n",
" elif fun(b) * fun(c) < 0:\n",
" a = c\n",
" msg = \"Maximum number of iterations reached.\"\n",
" return c, msg"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-231.6624790355399"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"raiz, _ = bisection(fun, -300.0, -200.0)\n",
"raiz"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-7.723180237198344e-15"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"raiz, _ = bisection(fun, -100.0, 0.100)\n",
"raiz"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"431.66247903554"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"raiz, _ = bisection(fun, 300.0, 500.0)\n",
"raiz"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"**Pregunta**\n",
"\n",
"Para al menos una de las raíces pruebe con diferentes rangos de localización de la raíz e indique el desempeño del algoritmo en términos del número de iteraciones.\n",
"\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Paso 2b: Afinación de la raíz por el método de Newton-Raphson"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recordemos el método de Newton-Raphson el cual consiste en extender la recta tangente en el punto actual $x_i$ hasta que cruce el cero para luego hacer la aproximación en la abscisa de dicho punto de cruce (ver figura). Considerando la pendiente en $x_i$ se tiene que:\n",
"\n",
"$$\n",
"f'(x_i)\\equiv m=\\frac{0-f(x_i)}{x_{i+1}-x_i}\n",
"$$\n",
"\n",
"resultando en la iteración de Newton-Raphson:\n",
"\n",
"$$\n",
"x_{i+1}=x_i-\\frac{f(x_i)}{f'(x_i)}\n",
"$$\n",
"\n",
"\n",
"
\n",
" \n",
"
\n",
"\n",
"En el siguiente bloque de código se implementa la iteración de Newton-Raphson. Note que además de la función, el algoritmo requiere como parámetros de entrada la primera derivada y un valor inicial de la raíz. El algoritmo puede fallar por alcanzar el máximo número de iteraciones permitido o por llegar a puntos en donde la derivada es cercana a cero."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Comentar apropiadamente el siguiente código\n",
"\n",
"def newton(fun, grad, x, niter=50, ftol=1e-8, verbose=False):\n",
" \"\"\"Use Newton method to estimate the root of a real function\n",
"\n",
" Parameters\n",
" ----------\n",
" fun : function\n",
" Function.\n",
" grad : function\n",
" Derivative of the function.\n",
" x : float\n",
" Initial estimate.\n",
" niter : int, optional\n",
" Maximum number of iterations.\n",
" ftol : float, optional\n",
" Tolerance for the root.\n",
" verbose : bool, optional\n",
" If True, prints each iteration.\n",
"\n",
" Returns\n",
" -------\n",
" x : array\n",
" Approximated root.\n",
" msg : string\n",
" Message regarding the success of the method.\n",
" \"\"\"\n",
" msg = \"Maximum number of iterations reached.\"\n",
" for cont in range(niter):\n",
" if abs(grad(x)) < ftol:\n",
" x = None\n",
" msg = \"Derivative near to zero.\"\n",
" break\n",
" if verbose:\n",
" print(\"n: {}, x: {}\".format(cont, x))\n",
" x = x - fun(x)/grad(x)\n",
" if abs(fun(x)) < ftol:\n",
" msg = \"Root found with desired accuracy.\"\n",
" break\n",
" return x, msg"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Probemos con diferentes puntos iniciales."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n: 0, x: -230.0\n",
"n: 1, x: -231.67883211678833\n",
"n: 2, x: -231.6624805928219\n"
]
},
{
"data": {
"text/plain": [
"(-231.66247903554, 'Root found with desired accuracy.')"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"newton(fun, deriv, -230.0, verbose=True)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n: 0, x: -51.0\n",
"n: 1, x: 10.94059640374946\n",
"n: 2, x: 0.20496836097227877\n",
"n: 3, x: 8.378324920088609e-05\n",
"n: 4, x: 1.4039249228220071e-11\n"
]
},
{
"data": {
"text/plain": [
"(3.9420326066980225e-25, 'Root found with desired accuracy.')"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"newton(fun, deriv, -51.0, verbose=True)\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(431.66247903554, 'Root found with desired accuracy.')"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"newton(fun, deriv, 430.0, verbose=False)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Glosario de términos"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**Raíz de una función $f(x)$:** Valor (o valores) de $x$ para los cuales se da la condición $f(x)=0$.\n",
"\n",
"**Convergencia:** Aproximaciones sucesivas a un valor buscado a través de un proceso iterativo.\n",
"\n",
"**Tolerancia:** Número real para definir el cero computacional en un algoritmo.\n",
"\n",
"**Método de bisección:** Algoritmo que localiza la raíz de una función $f(x)$ a través de particiones sucesivas del intervalo de búsqueda.\n",
"\n",
"**Método de Newton-Raphson:** Algoritmo que localiza la raíz de una función $f(x)$ a partir de la ecuación de la tangente en una aproximación inicial de la raíz.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Actividad para la clase"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"El método de Newton-Raphson resulta de la ecuación de la tangente en un punto $x_i $asumido como raíz de la función. Este método tiene como desventaja la necesidad de requerir la primera derivada de la función $f'(x)$, la cual no siempre esta disponible. Se puede formular también una iteración alternativa que evite el uso de $f'(x)$ si se hace uso de la recta secante entre 2 puntos $x_i$ y $x_{i+1}$ como se ilustra en la figura.\n",
"\n",
"
\n",
" \n",
"
\n",
"\n",
"Para derivar el método partimos de dos puntos $x_i$ y $x_{i-1}$, y encontramos la recta\n",
"secante que pasa por estos. Extendiendo la misma hasta su cruce con cero resulta en:\n",
"\n",
"$$m = \\frac{0 - f(x_i)}{x_{i + 1} - x_{i}}\\, ,$$\n",
"\n",
"y sabiendo que $m = (f(x_i) - f(x_{i - 1}))/(x_{i} - x_{i - 1})$, se tiene la siguiente iteración:\n",
"\n",
"$$x_{i + 1} = x_i - f(x_i) \\frac{x_{i} - x_{i - 1}}{f(x_i) - f(x_{i - 1})}\\, .$$\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"1. Modificar la rutina de Newton-Raphson de manera que determine la raíz de $f(x)$ usando el método de la secante.\n",
"\n",
"2. Para la función del ejemplo determinar sus raíces usando los métodos de Newton-Raphson y de la secante y compare su desempeño en términos de precisión y número de iteraciones.\n",
"\n",
"
\n",
"\n",
"\n",
"A continuación se brinda una plantilla para la actividad."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def secante(fun, x0, x1, niter=50, ftol=1e-8):\n",
" \"\"\"Use secant method to estimate the root of a real function\n",
"\n",
" Parameters\n",
" ----------\n",
" fun : function\n",
" Function.\n",
" x0 : float\n",
" Initial estimate.\n",
" x1 : float\n",
" Initial estimate.\n",
" niter : int, optional\n",
" Maximum number of iterations.\n",
" ftol : float, optional\n",
" Tolerance for the root.\n",
"\n",
" Returns\n",
" -------\n",
" x : array\n",
" Approximated root.\n",
" \"\"\"\n",
" for cont in range(niter): # Modificar\n",
" grad = 1 # Aproximacion al gradiente\n",
" x = x - fun(x)/grad\n",
" if abs(fun(x)) < ftol:\n",
" break\n",
" return x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Referencias"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Gómez, Juan., Sierra, César., Vergara, Juan., Sáenz, Mario., and Guarín-Zapata, Nicolás. (2018) [Notas de clase: Mecánica del medio continuo](https://github.com/AppliedMechanics-EAFIT/Notas-MMC/raw/master/notas_de_clase/notas_medios.pdf). Universidad EAFIT.\n",
"\n",
"2. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (1996). Numerical recipes in Fortran 90. Cambridge: Cambridge university press."
]
},
{
"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": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 16,
"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.7.3"
},
"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
}
![\"Esquema](\"img/biseccion.svg\"\n",)
\n",
"![\"Esquema](\"img/newton_iter.svg\"\n",)
\n",
"![\"Esquema](\"img/secante_iter.svg\"\n",)
\n",
"