Resolución numérica de ecuaciones¶

A continuación se resuelven detalladamente once problemas de resolución numérica de ecuaciones, provenientes del libro Numerical Methods in Engineering with Python de Jaan Kiusalaas (Problem Set 4.1).

Esta libreta interactiva hecha en Jupyter fue creada por Iván Sánchez y Samuel Ochoa.

Ejercicio 1¶

Usar el método de Newton-Raphson y una calculadora de cuatro funciones (solo operaciones de suma, resta, multiplicación y división) para calcular $\sqrt[3]{75}$ con precisión de cuatro cifras significativas.

Partimos sabiendo que el valor que buscamos es $x = \sqrt[3]{75}$, y realizamos un proceso de "despeje inverso", de manera que quede una expresión $f(x) = 0$ que se pueda expresar usando las cuatro operaciones básicas solamente:

$x = \sqrt[3]{75} ~~ \Rightarrow ~~ x^3 - 75 = 0 ~~ \Rightarrow ~~ x \cdot x \cdot x - 75 = 0$

Luego, necesitamos la expresión de la derivada de $f(x)$:

$f'(x) = \frac{d}{dx}(x^3 - 75) = 3x^2 ~~ \Rightarrow ~~ f'(x) = 3 \cdot x \cdot x$

Para aplicar el método de Newton-Raphson necesitamos una estimación inicial de x. Buscamos un número que elevado al cubo sea cercano a 75. Como 4³ = 64 y 5³ = 125, se toma $x_0 = 4$ por ser el entero cuyo cubo está más cercano a 75. Alternativamente se puede tomar un valor no entero como 4.2, por ejemplo.

A partir de esta estimación inicial, el método de Newton-Raphson consiste en generar nuevos valores de x hasta que se llegue a un límite de repeticiones o hasta que la variación de una estimación a otra esté por debajo de una tolerancia determinada. Cada iteración es generada con la fórmula:

$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

Como se quieren cuatro cifras significativas, se establece una tolerancia de 0.0001 y se aplica redondeo en el algoritmo. El número máximo de repeticiones se establece en un número arbitrario suficientemente grande, en este caso 1000, el valor predeterminado. A continuación se presenta la implementación de este procedimiento en Python:

Si se observa la gráfica de la función $f(x)$ se evidencia que la raíz de la función se encuentra cerca del valor aproximado que obtuvimos:

Ejercicio 2¶

Encontrar la raíz real positiva más pequeña de $x^3-3.23x^2-5.54x+9.84=0$ mediante el método de bisección.

Para hacer esto se requiere un intervalo inicial donde esté la raíz, una tolerancia que estableceremos en 0.0001 y un número máximo de repeticiones, en este caso como el ejercicio anterior se dejará en 1000 repeticiones.

Para el intervalo tomaremos como límite izquierdo el 0, porque el enunciado dice que se buscan números reales positivos; para el límite derecho analizaremos los puntos críticos de la función.

Para hallar los puntos críticos necesitamos la expresión de la derivada de $f(x)$:

$f'(x)= \frac{d}{dx}\left(x^3-3.23x^2-5.54x+9.84\right) = 3x^2-6.46x-5.54$

Los puntos críticos son aquellos en los que la derivada de la función es igual a cero:

$f'(x) = 0 \Rightarrow 3x^2-6.46x-5.54 = 0 \\ x = \frac{6.46 \pm \sqrt{(-6.46)^2 - 4(3)(-5.54)}}{2\cdot 300} \Rightarrow x_1 \approx -0.65708 \lor x_2 \approx 2.81041$

Como $x_1$ es negativo, se toma un número cercano a $x_2$ por debajo como límite superior. Entonces, nuestro intervalo de búsqueda será (0, 2.8)

Ya con todos los datos iniciales, a continuación se presenta la implementación de este procedimiento en Python y el resultado de su ejecución:

Si se observa la gráfica de la función $f(x)$ se evidencia que la raíz de la función se encuentra cerca del valor aproximado que obtuvimos:

Ejercicio 3¶

La raíz no nula, positiva, más pequeña de $\cosh(x)\cos(x) - 1 = 0$ se encuentra en el intervalo (4, 5). Calcular esta raíz mediante el método de Ridder.

Para empezar, se explica en qué consiste el método de Ridder. El método de Ridder es un algoritmo de búsqueda de raíces basado en el método de falsa posición y en el uso de una función exponencial para aproximar sucesivamente la raíz de una función contínua $f(x)$.

En resumen, se requiere que la raíz de $f(x)$ esté dentro de un intervalo inicial $(x_1,x_2)$ tal que $f(x_1)\cdot f(x_2)<0$. El método comienza evaluando la función en su punto medio $x_3 = \frac{x_1 + x_2}{2}$. La fórmula iterativa para generar la siguiente aproximación de la raíz es la siguiente (en el libro se explica cómo se obtiene):

$x_4 = x_3 \pm (x_3 - x_1)\cdot\frac{f(x_3)}{\sqrt{(f(x_3))^2 - f(x_1)\cdot f(x_2)}}$

Se usa el signo más (+) si $f(x_1) - f(x_2) > 0$, y el signo menos (-) si $f(x_1) - f(x_2) < 0$. Después de calcular $x_4$, se determina un nuevo intervalo para la raíz y se aplica el proceso de nuevo hasta que la diferencia entre dos valores sucesivos de $x_4$ sea insignificante (menor a una tolerancia preestablecida).

El método de Ridder tiene la propiedad de que si $x_1$ y $x_2$ encierran la raíz, entonces $x_4$ siempre estará dentro del intervalo $(x_1, x_2)$, haciendo al método muy confiable. Además, se puede demostrar que este método converge cuadráticamente, haciéndolo mucho más rápido que los métodos de la secante y de falsa posición, y no necesita calcular la derivada de la función.

El libro provee el código de una implementación del método de Ridder, el cual es mostrado a continuación:

El ejercicio pide encontrar la raíz positiva más pequeña de $\cosh(x)\cdot\cos(x) - 1 = 0$ en el intervalo (4,5). A continuación se traduce y modifica el código proveído por el libro para mostrar el proceso de aproximación aplicado a este ejercicio y permitir un mayor número de iteraciones:

Si se observa la gráfica de la función $f(x)$ se evidencia que la raíz de la función se encuentra cerca del valor aproximado que obtuvimos:

Ejercicio 4¶

Resolver el ejercicio 3 por el método de Newton-Raphson.

Se nos indica encontrar la raíz más pequeña de la función $\cosh(x)\cos(x) - 1$ en el intervalo (4, 5), esta vez usando el método de Newton-Raphson.

Idealmente se haría un estudio de los puntos críticos de la función usando su derivada, pero como se puede observar esta consiste de funciones trigonométricas e hiperbólicas que imposibilitan la resolución de la ecuación por métodos comunes:

$f'(x) = 0 \Rightarrow \sinh(x)\cos(x) - \cosh(x)\sin(x) = 0 \\ \Rightarrow \frac{\sinh(x)}{\cosh(x)}=\frac{\sin(x)}{\cos(x)} \Rightarrow \tanh(x)=\tan(x)$

Como no tenemos manera de verificar los puntos críticos matemáticamente, aprovechamos la gráfica de la función hecha en el ejercicio anterior. Como se puede ver, la función tiene un mínimo local muy cercano a $x = 4$, lo cual puede traer problemas (y si trae problemas, como se explica en el ejercicio 8). Por estos motivos usaremos un valor inicial de 5.

A continuación se presenta el código en Python para realizar lo que se pide. Se usan los valores predeterminados para la tolerancia ($1.0 \times 10^{-9}$) y el número de repeticiones (1000).

Este valor es igual al obtenido mediante el método de Ridder en el ejercicio anterior. De igual manera, si se observa la gráfica de la función $f(x)$ se evidencia que la raíz de la función se encuentra cerca del valor aproximado que obtuvimos:

Ejercicio 6¶

Determinar las dos raíces de $\sin(x) + 3\cos(x) - 2 = 0$ que se encuentran en el intervalo (-2, 2). Usar el método de Newton-Raphson.

Para buscar las dos raíces de la función $f(x)=\sin(x) + 3\cos(x) - 2$ usando el método de Newton-Raphson, se necesitan dos valores iniciales. Primero analizamos los puntos críticos de la función:

$f'(x) = 0 \Rightarrow \frac{d}{dx}\left(\sin(x) + 3\cos(x) - 2\right) = 0 \Rightarrow \cos(x)-3\sin(x) = 0 \\ \Rightarrow \cos(x)=3\sin(x) \Rightarrow \frac{1}{3} = \frac{\sin(x)}{\cos(x)} \Rightarrow \tan(x) = \frac{1}{3} \\ \Rightarrow x = \arctan\left(\frac{1}{3}\right) \approx 0.32175 $

Entonces, necesitamos un valor inicial a la izquierda del punto crítico y el otro a la derecha. Los extremos del intervalo servirán para esto.

A continuación se muestra la implementación en Python del método de Newton-Raphson aplicado con valores iniciales de -2 y 2 respectivamente, usando los valores predeterminados de la tolerancia y el número máximo de repeticiones:

Si se observa la gráfica de la función $f(x)$ se evidencia que las raíces de la función en el intervalo dado se encuentran cerca de los valores aproximados que obtuvimos:

Ejercicio 8¶

Dibujar una gráfica de $f(x) = \cosh(x) \cos(x) - 1$ en el rango $0 \leq x \leq 10$.

1. Verificar a partir de la gráfica que la raíz no nula, positiva, más pequeña de f(x) = 0 se encuentra en el intervalo (4, 5).

Debido a que la escala de la gráfica con el rango completo no permite apreciar bien el intervalo (4,5), se traza otra gráfica (abajo) con el intervalo $0 \leq x \leq 5$. Así se puede apreciar que en este intervalo hay dos raíces, $x = 0$ y $x \approx 4.73$.

La raíz en $x \approx 4.73$ se obtuvo en los ejercicios 3 y 4, mientras que se puede comprobar que $x = 0$ es una raíz y no una asíntota reemplazando en la expresión original: $\cosh(0)\cos(0) - 1 = 0 ~~ \Rightarrow 1 \cdot 1 - 1 = 0 ~~ \Rightarrow 0 = 0$

Como se pide la raíz no nula, positiva, más pequeña de la función, se tiene que $x \approx 4.73$ es la raíz que se busca, y pertenece al intervalo (4,5), como se aprecia en la gráfica.

1. Mostrar gráficamente que la fórmula de Newton-Raphson no convergería a esta raíz si se comienza con $x = 4$.

Como se observa en la gráfica de abajo, la recta tangente a $f'(x)$ en $x_0 = 4$ es bastante aplanada, debido a que está cerca de un mínimo local. La intersección de esta tangente con el eje horizontal, que sería el punto $x_1$ está muy alejada de $x_0$, en $x \approx 10.6629$.

Si se sigue iterando el proceso de calcular $x_i$ con la fórmula de Newton-Raphson, se obtiene $x \approx 10.9956$, una raíz de la función mucho mayor a la esperada $x \approx 4.73$.

Entonces, se confirma que si se inicia en $x = 4$ no se llega a la raíz no nula positiva más pequeña, debido a la cercanía de dicho valor inicial con un mínimo de la función.

Ejercicio 9¶

La ecuación $x^3 - 1.2x^2 - 8.19x + 13.23 = 0$ tiene una doble raíz cerca de $x = 2$. Determinar esta raíz con el método Newton-Raphson dentro de cuatro espacios decimales.

"Doble raíz" significa que la raíz en el punto que buscamos es también un punto crítico. Los puntos críticos de la función se obtienen igualando su derivada a cero:

$f'(x) = 0 \Rightarrow \frac{d}{dx}\left(x^3 - 1.2x^2 - 8.19x + 13.23\right) = 0 \Rightarrow 3x^2-2.4x-8.19 = 0 \\ x = \frac{2.4 \pm \sqrt{(-2.4)^2-4(3)(-8.19)}}{2\cdot 3} \Rightarrow x_1 = -1.3 \lor x_2 = 2.1 $

El punto crítico $x_2$ es el más cercano a $x = 2$, por lo que tenemos que comprobar que el método de Newton-Raphson produzca un valor aproximadamente igual a él. Como valor inicial tomamos un punto intermedio a los puntos críticos, como $x = 0$.

A continuación se implementa en Python el método de Newton-Raphson modificado para trabajar con cuatro decimales de precisión:

Así se comprueba que la raíz de la función se encuentra en $x = 2.1$, y es una doble raíz porque coincide con un punto crítico. A continuación se muestra la gráfica de $f(x)$, donde se aprecia que la raíz en $x = 2.1$ es a su vez un mínimo local:

Ejercicio 10¶

Escribir un programa que calcule todas las raíces de $f(x) = 0$ en un intervalo dado con el método de Ridder. Utilizar las funciones rootsearch y ridder. Puede usar el programa en el Ejemplo 4.3 como modelo. Probar el programa encontrando las raíces de $x \cdot \sin( x ) + 3\cos( x ) - x = 0$ en (-6, 6).

El libro provee una función llamada rootsearch, la cual busca una raíz de la función f(x) en el intervalo (a,b) en incrementos de dx. Su resultado es un tuple (x1,x2) cuyos elementos son los extremos del intervalo donde se encuentra la raíz, si la búsqueda es exitosa. Si en lugar de números, los valores de x1 y x2 son None, significa que no se encontró ninguna raíz.

Después de que la primera raíz (la más cercana al extremo izquierdo) ha sido detectada, rootsearch puede ser llamada de nuevo con x2 como extremo izquierdo para encontrar la siguiente raíz. Esto puede ser repetido siempre que la función detecte una raíz.

El código de rootsearch, tal como aparece en el libro, es el siguiente:

Para encontrar todas las raíces de $f(x) = x \cdot \sin( x ) + 3\cos( x ) - x = 0$ en (-6,6) usando el método de Ridder podemos aprovechar esta función rootsearch para generar los intervalos que verificará el método de Ridder. A continuación se implementa este procedimiento en Python:

Si se observa la gráfica de la función $f(x)$ se evidencia que la función tiene raíces aproximadamente en los puntos encontrados por el algoritmo:

Ejercicio 11¶

Repetir el ejercicio 10 con el método de Newton-Raphson.

Se pide crear un programa que encuentre todas las raíces de $f(x) = x \cdot \sin( x ) + 3\cos( x ) - x = 0$ en (-6, 6), esta vez usando el método de Newton-Raphson.

Primero que nada se buscan los puntos críticos de la función:

$f'(x) = 0 \Rightarrow \frac{d}{dx}\left(x \cdot \sin( x ) + 3\cos( x ) - x\right) = 0 \\ \Rightarrow \sin(x) + x\cdot\cos(x) - 3\sin(x) - 1 = 0 \\ \Rightarrow x\cdot\cos(x) - 2\sin(x) - 1 = 0 $

No es evidente una forma de resolver esta ecuación algebraicamente, así que nos basaremos en la gráfica de $f(x)$ realizada en el ejercicio anterior. A partir de ahí vemos tres puntos críticos, un mínimo cerca de $x = -4$, un máximo cerca de $x = -1$, y otro mínimo cerca de $x = 4$; estos son valores que tenemos que evitar para que el método de Newton-Raphson pueda converger.

Estos puntos críticos (aproximados) dividen el intervalo (-6, 6) en cuatro partes: una entre -6 y -4, una entre -4 y -1, una entre -1 y 4, y una entre 4 y 6. Tomaremos un valor intermedio de cada subintervalo como valor inicial para el procedimiento, por ejemplo: $x = -5$, $x = -2$, $x = 0$, $x = 5$.

A continuación se implementa en Python el método de Newton-Raphson y se aplica a estos valores iniciales, utilizando los valores predeterminados para la tolerancia y el número máximo de repeticiones:

Las tres raíces obtenidas que pertenecen al intervalo (-6, 6) son iguales a las obtenidas en el ejercicio anterior. De igual manera, si se observa la gráfica de la función $f(x)$ se evidencia que las raíces de la función se encuentran cerca de los valores aproximados que obtuvimos:

Ejercicio 13¶

Calcular todas las raíces reales positivas de $x^4 + 2x^3 - 7x^2 + 3 = 0$

Para hallar las raíces de $f(x) = x^4 + 2x^3 - 7x^2 + 3 = 0$ podemos usar el método de Newton-Raphson, pues no requiere encerrar las raíces en intervalos, sino dar estimaciones iniciales.

Las estimaciones iniciales no pueden ser puntos críticos, por lo que tenemos que identificar estos puntos igualando la derivada de la función a cero:

$f'(x) = 0 \Rightarrow \frac{d}{dx}\left(x^4 + 2x^3 - 7x^2 + 3\right) = 0 \\ \Rightarrow 4x^3 + 6x^2 - 14x = 0 \Rightarrow 2x(2x^2+3x-7) = 0 \\ \Rightarrow x_1 = 0 ~~ \lor ~~ 2x^2 + 3x - 7 = 0 \\ \Rightarrow x = \frac{-3 \pm \sqrt{3^2-4(2)(-7)}}{2\cdot 2} \Rightarrow x_2 \approx -2.76556 \lor x_3 \approx 1.26556$

La existencia de tres puntos críticos significa que se puede dividir el dominio en cuatro intervalos: entre $-\infty$ y -2.76556, entre -2.76556 y 0, entre 0 y 1.26556, y entre 1.26556 y $\infty$. Vamos a tomar un punto de cada intervalo como valores iniciales, en este caso: $x = -3$, $x = -1$, $x = 1$ y $x = 2$.

A continuación se muestra la implementación del método de Newton-Raphson que venimos manejando, aplicada a cada uno de estos cuatro valores iniciales, con la tolerancia y el número máximo de repeticiones predeterminados.

Si se observa la gráfica de la función $f(x)$ se evidencia que las raíces de la función se encuentran cerca de estos cuatro valores aproximados que obtuvimos:

Ejercicio 18¶

La ecuación de Bernoulli para el flujo de fluido en un canal abierto con una pequeña protuberancia es: $\frac{Q^2}{ 2gb^2{h_0}^2 } + h_0 = \frac{Q^2}{ 2gb^2 h^2 } + h + H$ donde:

• $Q = 1.2 m^3/s$ es la tasa de flujo volumétrico.
• $g = 9.81 m/s^2$ es la aceleración gravitatoria.
• $b = 1.8 m$ es la anchura del canal.
• $h_0 = 0.6 m$ es el nivel de agua aguas arriba.
• $H = 0.075 m$ es la altura de la protuberancia.
• $h = ?$ es el nivel del agua sobre la protuberancia.

Determinar h.

Para hallar el valor de h se puede pasar todo a un solo lado de la ecuación, de forma que quede una expresión $f(h) = 0$ a la cual se le puede aplicar el método de Newton-Raphson:

$\frac{Q^2}{ 2gb^2{h_0}^2 } + h_0 = \frac{Q^2}{ 2gb^2 h^2 } + h + H \\ \frac{Q^2}{ 2gb^2 h^2 } + h + H - \frac{Q^2}{ 2gb^2{h_0}^2 } - h_0 = 0 \\ \Rightarrow f(h) = 0 $

Como esta función posee una expresión racional con $h^2$ en el denominador, posee una asíntota en $h = 0$. Al tratarse de la altura de una protuberancia, h no puede tomar valores negativos, así que se busca una raíz real no nula y positiva.

A continuación se analizan los puntos críticos de la función para determinar cuántos valores iniciales se tienen que tomar. Los puntos críticos se obtienen resolviendo la ecuación siguiente:

$f'(h) = 0 \Rightarrow \frac{d}{dx}\left(\frac{Q^2}{ 2gb^2}\cdot h^{-2} + h + H - \frac{Q^2}{ 2gb^2{h_0}^2 } - h_0\right) = 0 \\ \Rightarrow -2\cdot\frac{Q^2}{ 2gb^2 }\cdot h^-3 + 1 = 0 \Rightarrow 1-\frac{Q^2}{ gb^2h^3 } = 0 \\ \Rightarrow 1 = \frac{Q^2}{ gb^2h^3 } \Rightarrow h^3 = \frac{Q^2}{ gb^2 } \Rightarrow h = \sqrt[3]{\frac{Q^2}{ gb^2 }} = \sqrt[3]{\frac{(1.2 m^3/s)^2}{(9.81 m/s^2)(1.8m)^2}} \approx 0.35649 m $

Este valor de h se debe evitar, así que tomamos un valor inicial entre 0 y 0.35649, y otro entre 0.35649 y $\infty$. En este caso tomaremos $x = 0.2$ y $x = 0.6$.

A continuación se implementa en Python el método de Newton-Raphson para $f(x)$ y se aplica con los valores iniciales mencionados, y usando los valores predeterminados para la tolerancia y el número máximo de iteraciones:

Si se observa la gráfica de la función $f(h)$ se evidencia que las raíces de la función se encuentran cerca de estos dos valores aproximados que obtuvimos:

La primera raíz encontrada parece estar alejada de la posición real de la raíz en la gráfica, lo cual indica que hay una pérdida de precisión, posiblemente causada porque las fórmulas tienen varias divisiones y restas, que se acentúa conforme los valores con los que se trabaja son más pequeños.