🔍 Resolución de Ecuaciones No Lineales con Python #
¡Bienvenido a la primera cápsula de la serie de Métodos Numéricos! En esta sección, nos adentraremos en uno de los desafíos más comunes en la ingeniería: encontrar las raíces de ecuaciones no lineales. A menudo, las ecuaciones que describen sistemas físicos no tienen una solución analítica sencilla, obligándonos a recurrir a algoritmos numéricos para encontrar una solución aproximada.
En esta cápsula, te guiaremos a través de los conceptos teóricos, crearemos nuestros propios algoritmos desde cero, y finalmente, te mostraremos cómo las librerías de Python ya tienen herramientas listas para usar, haciéndote más eficiente.
Conceptos Teóricos: El Problema de la Raíz #
Una raíz de una función \(f(x)\) es un valor \(x\) tal que \(f(x) = 0\). En términos gráficos, es el punto donde la función cruza el eje horizontal. En ingeniería, esto puede representar, por ejemplo, el punto de equilibrio de un sistema, el instante en que la velocidad de un objeto se anula, o la frecuencia de resonancia de una estructura.
La búsqueda de estas raíces se basa en métodos iterativos, donde se comienza con una estimación inicial y se refina progresivamente hasta que se cumple un criterio de convergencia, como un error menor a una tolerancia definida.
Método de Punto Fijo (Iteración Simple) #
El método de Punto Fijo es la forma más intuitiva de un proceso iterativo y la base de muchos otros métodos. En lugar de resolver \(f(x) = 0\) directamente, el método requiere que la función se reescriba en una forma equivalente:
$$x = g(x)$$Un valor \(p\) es un punto fijo de \(g(x)\) si \(p = g(p)\). Si \(p\) es un punto fijo de \(g(x)\), entonces también es una raíz de \(f(x) = x - g(x) = 0\).
El algoritmo es simple: se toma una estimación inicial \(x_0\) y se genera una secuencia de aproximaciones usando la relación de recurrencia:
$$x_{n+1} = g(x_n)$$Criterio de Convergencia #
La convergencia de este método no siempre está garantizada. Una condición suficiente para que la iteración converja a un punto fijo único \(p\) en un intervalo \([a, b]\) es que:
- \(g(x)\) mapee el intervalo \([a, b]\) en sí mismo (es decir, \(g(x) \in [a, b]\) para todo \(x \in [a, b]\)).
- La magnitud de la derivada sea menor que 1 en ese intervalo: \(|g'(x)| < 1\) para todo \(x \in (a, b)\).
Si la segunda condición se cumple, se dice que la función \(g(x)\) es una contracción en el intervalo, garantizando la convergencia.
Pseudocódigo: Método de Punto Fijo #
funcion punto_fijo(g, x0, tol, max_iter):
x_n = x0
PARA i DESDE 1 HASTA max_iter HACER
x_next = g(x_n)
SI abs(x_next - x_n) < tol ENTONCES
RETORNAR x_next
FIN SI
x_n = x_next
FIN PARA
RETORNAR "No converge"
Diagrama de Bloques: Método de Punto Fijo #
graph TD;
A[Inicio] --> B["Establecer x_n = x0"];
B --> C{"Bucle: Para i = 1 hasta max_iter"};
C --> D["Calcular x_next = g(x_n)"];
D --> E{"¿abs(x_next - x_n) < tol?"};
E -- Sí --> F["Raíz encontrada: x_next"];
E -- No --> I["Actualizar x_n = x_next"];
I --> C;
F --> J[Fin];
Método de Bisección #
El método de Bisección es uno de los más simples y robustos para encontrar raíces. Se basa en el Teorema del Valor Intermedio, que establece que si una función continua \(f(x)\) tiene signos opuestos en los extremos de un intervalo \([a, b]\), es decir, \(f(a)f(b) < 0\), entonces debe existir al menos una raíz en ese intervalo.
El algoritmo funciona de la siguiente manera:
- Se elige un intervalo inicial \([a, b]\) donde se sabe que hay una raíz.
- Se calcula el punto medio \(c = (a+b)/2\).
- Se evalúa \(f(c)\). Si \(f(c) = 0\), se encontró la raíz.
- Si no, se decide si la raíz está en \([a, c]\) o en \([c, b]\) analizando el signo de \(f(a)f(c)\). Si es negativo, la raíz está en \([a, c]\); de lo contrario, está en \([c, b]\).
- Se repite el proceso con el nuevo intervalo, el cual es la mitad de grande que el anterior, hasta que la longitud del intervalo sea menor que la tolerancia deseada.
Pseudocódigo: Método de Bisección #
Se dan los valores a y b
REPITE INDEFINIDAMENTE
c ← (a + b) / 2
SI abs(f(c)) < tol ENTONCES
RETORNAR c
FIN SI
SI signo(f(a)) != signo(f(c)) ENTONCES
b ← c
SINO
a ← c
FIN SI
FIN REPETIR
Diagrama de Bloques: Método de Bisección #
graph TD;
A[Inicio] --> B["Elegir intervalo [a, b] tal que f(a)f(b) < 0"];
B --> C{Iteración};
C --> D["Calcular punto medio c = (a + b) / 2"];
D --> E{"¿f(c) < tol?"};
E -- Sí --> F["Raíz encontrada: c"];
E -- No --> J{"f(a)f(c) < 0?"};
J -- Sí --> K["Nuevo intervalo: [a, c]"];
J -- No --> L["Nuevo intervalo: [c, b]"];
K ---> C
L ---> C
F --> O[Fin];
Método de Newton-Raphson #
El método de Newton-Raphson es mucho más rápido y eficiente que el de Bisección, pero requiere que la función sea diferenciable y que se conozca su derivada. Se basa en la idea de usar la tangente de la función en un punto para aproximar la siguiente iteración.
La fórmula de recurrencia es:
$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$Donde:
- \(x_n\) es la aproximación actual.
- \(f(x_n)\) es la función evaluada en \(x_n\).
- \(f'(x_n)\) es la derivada de la función evaluada en \(x_n\).
Este método converge mucho más rápido que el de Bisección, pero puede fallar si la derivada es cero o muy pequeña, o si la estimación inicial está muy lejos de la raíz.
Pseudocódigo: Método de Newton-Raphson #
x_n = x0
PARA i DESDE 1 HASTA max_iter HACER
f_val = f(x_n)
df_val = df(x_n)
x_next = x_n - f_val / df_val
SI abs(x_next - x_n) < tol ENTONCES
RETORNAR x_next
FIN SI
x_n = x_next
FIN PARA
Diagrama de Bloques: Método de Newton-Raphson #
graph TD;
A[Inicio] --> B["Establecer x_n = x0"];
B --> C{"Bucle: Para i = 1 hasta max_iter"};
C --> D["Calcular f_val = f(x_n)"];
D --> E["Calcular df_val = df(x_n)"];
E --> F["Calcular x_next = x_n - f_val / df_val"];
F --> G{"¿abs(x_next - x_n) < tol?"};
G -- Sí --> H["Retornar x_next"];
G -- No --> I["Actualizar x_n = x_next"];
I --> C;
H --> J[Fin];
Implementando los Métodos con Python #
A continuación, crearemos las funciones en Python para resolver ecuaciones utilizando los métodos de Bisección y Newton-Raphson. Utilizaremos NumPy para los cálculos numéricos y SymPy para encontrar la derivada de la función de forma simbólica, combinando las habilidades de las cápsulas anteriores.
import numpy as np
import sympy as sp
# ----------------- Funciones para los Métodos -----------------
def punto_fijo(g, x0, tol=1e-6, max_iter=100):
"""
Resuelve una ecuación no lineal usando el método de Punto Fijo.
Args:
g (function): La función de iteración g(x) = x.
x0 (float): Estimación inicial.
tol (float): Tolerancia del error.
max_iter (int): Número máximo de iteraciones.
Returns:
float: La raíz aproximada.
"""
x_n = x0
for _ in range(max_iter):
x_next = g(x_n)
if np.abs(x_next - x_n) < tol:
return x_next
x_n = x_next
print("El método de Punto Fijo no convergió en el número máximo de iteraciones.")
return None
def biseccion(f, a, b, tol=1e-6, max_iter=100):
"""
Resuelve una ecuación no lineal usando el método de Bisección.
Args:
f (function): La función a la cual se le busca la raíz.
a (float): Extremo inferior del intervalo.
b (float): Extremo superior del intervalo.
tol (float): Tolerancia del error.
max_iter (int): Número máximo de iteraciones.
Returns:
float: La raíz aproximada.
"""
if f(a) * f(b) >= 0:
print("El método de Bisección no garantiza una raíz en este intervalo.")
return None
iter_count = 0
while (b - a) / 2 > tol and iter_count < max_iter:
c = (a + b) / 2
if f(c) == 0:
return c
elif f(a) * f(c) < 0:
b = c
else:
a = c
iter_count += 1
return (a + b) / 2
def newton_raphson(f, x0, tol=1e-6, max_iter=100):
"""
Resuelve una ecuación no lineal usando el método de Newton-Raphson.
Args:
f (function): La función a la cual se le busca la raíz.
x0 (float): Estimación inicial de la raíz.
tol (float): Tolerancia del error.
max_iter (int): Número máximo de iteraciones.
Returns:
float: La raíz aproximada.
"""
x = sp.symbols('x')
f_sym = f(x)
f_prime_sym = sp.diff(f_sym, x)
f_numeric = sp.lambdify(x, f_sym, 'numpy')
f_prime_numeric = sp.lambdify(x, f_prime_sym, 'numpy')
x_n = x0
for _ in range(max_iter):
f_val = f_numeric(x_n)
f_prime_val = f_prime_numeric(x_n)
if np.abs(f_prime_val) < 1e-9: # Evitar división por cero
print("Derivada muy pequeña. El método no converge.")
return None
x_next = x_n - f_val / f_prime_val
if np.abs(x_next - x_n) < tol:
return x_next
x_n = x_next
print("El método no convergió en el número máximo de iteraciones.")
return None
# ----------------- Ejemplos de Uso -----------------
# Definir la función f(x) = x^3 - 2x - 5
def mi_funcion(x):
return x**3 - 2*x - 5
# Para Punto Fijo, reescribimos f(x)=0 como x = g(x). Por ejemplo: x = (2x + 5)^(1/3)
def mi_funcion_g(x):
return (2*x + 5)**(1/3)
# Resolver usando Punto Fijo
raiz_pf = punto_fijo(mi_funcion_g, 2.0)
print(f"Raíz encontrada con Punto Fijo: {raiz_pf:.6f}")
# Resolver usando Bisección
raiz_biseccion = biseccion(mi_funcion, 2, 3)
print(f"Raíz encontrada con Bisección: {raiz_biseccion:.6f}")
# Resolver usando Newton-Raphson
raiz_newton = newton_raphson(mi_funcion, 2)
print(f"Raíz encontrada con Newton-Raphson: {raiz_newton:.6f}")
Funciones de Librería: Eficiencia al Máximo #
Si bien es importante entender y poder implementar los algoritmos, en la práctica profesional, los ingenieros utilizan librerías optimizadas para resolver sistemas lineales de manera rápida y eficiente. SciPy, que ya vimos en una cápsula anterior, tiene un módulo llamado optimize diseñado específicamente para este tipo de problemas.
import numpy as np
from scipy import optimize
# Definir la función para usar con SciPy
def mi_funcion_scipy(x):
return x**3 - 2*x - 5
# Método de Bisección de SciPy
raiz_biseccion_scipy = optimize.bisect(mi_funcion_scipy, 2, 3)
print(f"Raíz encontrada con SciPy bisect: {raiz_biseccion_scipy:.6f}")
# Método de Newton de SciPy
# Requiere una estimación inicial y la función derivada (opcional pero recomendada)
raiz_newton_scipy = optimize.newton(mi_funcion_scipy, 2)
print(f"Raíz encontrada con SciPy newton: {raiz_newton_scipy:.6f}")
Casos Aplicados en Ingeniería Mecánica #
1. Diseño de un Cilindro a Presión #
El diseño de recipientes a presión es crucial en la ingeniería. Para un cilindro de paredes gruesas sometido solamente a presión interna (\(P_i\)), los esfuerzos en cualquier punto \(r\) está dado por:
$$\sigma_1 = \sigma_t = \frac{r_i^2P_i}{r_o^2-r_i^2}\left(1 + \frac{r_o^2}{r^2}\right)$$$$ \sigma_2 = \sigma_r = \frac{r_i^2P_i}{r_o^2-r_i^2}\left(1 - \frac{r_o^2}{r^2}\right)$$$$ \sigma_3 = \sigma_l = \frac{P_ir_i^2}{r_o^2-r_i^2}$$
Con estos esfuerzos poodemos calcualer el esfuerzo equivalente de Von Mises.
$$\sigma_{v} = \sqrt{\frac{(\sigma_1 - \sigma_2)^2 - (\sigma_1 - \sigma_3)^2 - (\sigma_2 - \sigma_3)^2}{2}}$$Problema: Queremos encontrar el radio exterior (\(r_e\)) que resulte en un esfuerzo de Von Mises máximo (\(\sigma_{v,máx}\)) de \(200 \text{ MPa}\) en la pared interior (\(r=r_i\)).
Datos:
- Presión interna, \(P_i = 10 \text{ MPa}\)
- Radio interior, \(r_i = 0.5 \text{ m}\)
- Esfuerzo máximo deseado, \(\sigma_{v,máx} = 200 \text{ MPa}\)
Para resolver esto, primero debemos reorganizar la fórmula del esfuerzo de Von Mises para que sea una ecuación no lineal igual a cero, \(f(r_e) = 0\).
$$f(r_e) = \sigma_{v,máx} - \sigma_{v} = 0$$Ahora, implementamos esta función en Python y usamos el método de Newton-Raphson para encontrar \(r_e\).
from scipy import optimize
# Definir la función para resolver
def von_mises_func(re):
Pi = 10
ri = 0.5*1000
sigma_max = 200
re = re*1000 # Convertir a mm
sigma_1 = ri**2*Pi/(re**2 - ri**2) * (1 + re**2/re**2)
sigma_2 = ri**2*Pi/(re**2 - ri**2) * (1 - re**2/re**2)
sigma_3 = Pi*ri**2/(re**2 - ri**2)
# La ecuación es: f(re) = sigma_max - von_mises(re) = 0
return sigma_max - ((sigma_1-sigma_2)**2 + (sigma_2-sigma_3)**2 + (sigma_3-sigma_1)**2)**0.5 / 2**0.5
print(von_mises_func(0.52))
# Usar el método de Newton-Raphson de SciPy para encontrar la raíz
# Se necesita una estimación inicial. Probaremos con 0.501 m.
re_sol = optimize.newton(von_mises_func, 0.501)
print(f"El radio exterior requerido es: {re_sol:.4f} metros")
2. Análisis de un Circuito Hidráulico Simple #
En el análisis de fluidos, las ecuaciones a menudo involucran términos no lineales que modelan la pérdida de carga por fricción. La ecuación que describe la presión en un punto de un circuito hidráulico puede ser:
$$P = P_{fuente} - K Q^2 - \Delta P_{válvula} = 0$$Donde \(Q\) es el flujo volumétrico, \(K\) es una constante de fricción y \(\Delta P_{válvula}\) es la caída de presión en una válvula que depende de \(Q\). Supongamos que esta última es \(\Delta P_{válvula} = 50 Q + 100\).
Problema: En un sistema hidráulico, se tiene una presión de fuente de \(1500 \text{ kPa}\) y una constante de fricción \(K = 100 \text{ kPa} / (\text{m}^3/\text{s})^2\). Queremos encontrar el flujo volumétrico \(Q\) en el punto de equilibrio donde la presión final es cero.
La ecuación se puede reescribir como:
$$f(Q) = 1500 - 100 Q^2 - (50 Q + 100) = 0$$Simplificando, obtenemos una ecuación cuadrática que se podría resolver analíticamente, pero la usaremos para demostrar el uso de los métodos numéricos.
$$f(Q) = -100 Q^2 - 50 Q + 1400 = 0$$Usaremos el método de Bisección para encontrar la raíz positiva del flujo volumétrico.
from scipy import optimize
# Definir la función para el problema hidráulico
def circuito_func(Q):
return -100 * Q**2 - 50 * Q + 1400
# Usar el método de Bisección de SciPy
# Se sabe que la raíz positiva está en el intervalo [0, 5]
Q_sol = optimize.bisect(circuito_func, 0, 5)
print(f"El flujo volumétrico es: {Q_sol:.4f} m^3/s")
Ejercicios Prácticos #
Pon a prueba lo aprendido con los siguientes 20 ejercicios.
Parte 1: Implementación y Conceptos (Ejercicios 1-10) #
- Función: Encuentra la raíz de \(f(x) = x^2 - 4\) usando el método de Bisección en el intervalo \([1, 3]\).
- Función: Usa el método de Newton-Raphson para encontrar la raíz de \(f(x) = x^2 - 4\) con una estimación inicial de \(x_0 = 3\).
- Teoría: ¿Qué problema puede ocurrir con el método de Newton-Raphson si la derivada de la función es cercana a cero?
- Función: Encuentra la raíz de \(f(x) = \sin(x) - 0.5\) en el intervalo \([0, \pi/2]\) usando Bisección.
- Función: Usa Newton-Raphson para encontrar la raíz del ejercicio 4 con \(x_0 = 0.5\).
- Teoría: Explica por qué el método de Bisección siempre converge, a diferencia del método de Newton-Raphson.
- Función: Encuentra la raíz de \(f(x) = e^x - 3x\) en el intervalo \([1, 2]\) usando Bisección.
- Función: Usa Newton-Raphson para encontrar la raíz del ejercicio 7 con \(x_0 = 1.5\).
- Función: Modifica la función de Newton-Raphson para que también retorne el número de iteraciones que tomó.
- Función: Modifica la función de Bisección para que retorne una lista de todos los puntos medios calculados en el proceso.
Parte 2: Aplicaciones y Uso de Librerías (Ejercicios 11-20) #
- Aplicación: La ecuación para la deformación de una barra bajo una carga esférica es \(f(x) = \frac{x^3}{3} - 4\). Encuentra la deformación \(x\) para \(f(x) = 0\) usando
scipy.optimize.bisecten el intervalo \([2, 3]\). - Aplicación: El factor de fricción \(f\) en tuberías se puede encontrar con la ecuación de Colebrook. Para un caso particular, la ecuación es \(\frac{1}{\sqrt{f}} + 2 \log_{10}\left(\frac{2.51}{Re\sqrt{f}}\right) = 0\). Usa
scipy.optimize.newtonpara encontrar \(f\) si \(Re = 10^5\) y la estimación inicial es \(f_0 = 0.01\). - Aplicación: La temperatura de un reactor químico se modela con \(T(t) = 100(1 - e^{-0.1t}) + 20\). ¿En qué momento \(t\) la temperatura será de \(85 °C\)? Define la función \(f(t) = T(t) - 85\) y usa Bisección para encontrar la raíz en el intervalo \([1, 20]\).
- Aplicación: El ángulo \(\theta\) que un satélite debe tener para una trayectoria esférica es \(\frac{1}{\sin(\theta)} + \cos(\theta) = 3\). Encuentra \(\theta\) en el intervalo \([0.5, 1.5]\) usando
scipy.optimize.bisect. - Aplicación: La ecuación de estado de Van der Waals para un gas es \(P = \frac{RT}{v-b} - \frac{a}{v^2}\). Para valores dados de \(P, R, T, a, b\), encuentra el volumen molar \(v\). Define la función \(f(v) = P - (\frac{RT}{v-b} - \frac{a}{v^2})\) y usa Newton-Raphson para encontrar \(v\).
- Aplicación: Un problema de vibraciones mecánicas da la siguiente ecuación característica: \(f(\omega) = \omega^2 - \frac{k}{m}\). Encuentra la frecuencia natural \(\omega\) para \(k=100 N/m\) y \(m=5 kg\) usando Bisección.
- Aplicación: El diseño de una leva requiere encontrar la raíz de \(f(x) = 1.5 - x - \sin(x)\). Usa
scipy.optimize.newtonpara encontrar la raíz con una estimación inicial de \(x_0 = 1\). - Aplicación: Una ecuación de transferencia de calor es \(f(T) = hA(T - T_{amb}) - \epsilon\sigma A T^4\). Encuentra la temperatura de equilibrio \(T\) cuando el término radiativo (\(\epsilon\sigma A T^4\)) iguala a la transferencia por convección. Usa Newton-Raphson para encontrar la raíz para valores dados.
- Análisis: Compara la velocidad de convergencia de los métodos de Bisección y Newton-Raphson para el problema 13. ¿Cuál de los dos métodos converge más rápido?
- Reto: Resuelve el problema 12 usando tu propia función de Newton-Raphson. Para ello, necesitarás encontrar la derivada de la ecuación de Colebrook, ya sea simbólicamente con SymPy o manualmente.