🎯 Módulo 8: Newton-Raphson Generalizado
#

🧩 Unificando la Solución de Ecuaciones
#

El Método de Newton-Raphson Generalizado extiende el método de Newton para resolver sistemas de $n$ ecuaciones no lineales con $n$ incógnitas. Es el método iterativo más eficiente, combinando cálculo (el Jacobiano) con álgebra lineal (solución de sistemas).

Buscamos el vector solución $\mathbf{x} = (x_1, x_2, \dots, x_n)$ que satisfaga el sistema $\mathbf{F}(\mathbf{x}) = \mathbf{0}$.

🧠 Fundamento Teórico
#

La fórmula iterativa se basa en la aproximación de Taylor de primer orden:

$$\mathbf{x}^{k+1} = \mathbf{x}^k - \mathbf{J}(\mathbf{x}^k)^{-1} \mathbf{F}(\mathbf{x}^k)$$

En la práctica, la fórmula se resuelve como un Sistema de Ecuaciones Lineales para encontrar el vector de corrección $\boldsymbol{\Delta}\mathbf{x}^k$:

$$\mathbf{J}(\mathbf{x}^k) \cdot \boldsymbol{\Delta}\mathbf{x}^k = -\mathbf{F}(\mathbf{x}^k)$$

Donde $\mathbf{x}^{k+1} = \mathbf{x}^k + \boldsymbol{\Delta}\mathbf{x}^k$.

La Matriz Jacobiana ($\mathbf{J}$)
#

$\mathbf{J}$ es la matriz de derivadas parciales que contiene todas las tasas de cambio del sistema:

$$\mathbf{J}(\mathbf{x}) = \begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1} & \cdots & \frac{\partial f_n}{\partial x_n} \end{pmatrix}$$

💻 Ejemplo de Programación Manual con Auxiliares (Jacobiano Asistido)
#

Este ejemplo ilustra el algoritmo de Newton paso a paso, usando scipy.differentiate.jacobian para calcular la matriz $\mathbf{J}$ numéricamente y numpy.linalg.solve para resolver el sistema lineal, delegando las tareas pesadas a las librerías.

Sistema: $f_1(x, y) = x^2 + y^2 - 4 = 0$ y $f_2(x, y) = e^x - y - 1 = 0$.

import numpy as np
import math
# Función especializada para el cálculo del Jacobiano Numérico
from scipy.differentiate import jacobian 

# 1. Función Residual F(x)
def F_vector(x):
    """Calcula el vector de residuos F(x)."""
    f1 = x[0]**2 + x[1]**2 - 4
    f2 = np.exp(x[0]) - x[1] - 1
    return np.array([f1, f2])

def newton_generalizado_asistido(x_guess, max_iter=10, tol=1e-6):
    
    x_k = np.array(x_guess, dtype=float)
    
    print("--- Proceso Newton-Raphson Generalizado (Jacobiano Asistido) ---")
    
    for k in range(1, max_iter + 1):
        
        # 1. Calcular el Vector F (Residuos)
        F = F_vector(x_k)

        # 2. Calcular la Matriz Jacobiana J usando scipy.differentiate.jacobian
        J_result = jacobian(F_vector, x_k)
        J = J_result.df

        # 3. Resolver J * delta_x = -F usando NumPy
        menos_F = -F
        try:
            # np.linalg.solve: Solución del sistema lineal J * delta_x = -F
            delta_x = np.linalg.solve(J, menos_F)
        except np.linalg.LinAlgError:
            print(f"Error: La Matriz Jacobiana es singular en iteración {k}.")
            return x_k

        # 4. Criterio de Convergencia
        error = np.linalg.norm(delta_x)
        print(f"Iteración {k}: x={x_k}, Error (norma de corrección)={error:.8e}")
        
        if error < tol:
            print(f"Convergencia lograda en {k} iteraciones.")
            return x_k

        # 5. Actualización
        x_k += delta_x

    print("Máximo de iteraciones alcanzado. Sin convergencia.")
    return x_k

# Ejecución
initial_guess = [1.0, 1.0]
solution = newton_generalizado_asistido(initial_guess)
print(f"\nSolución final encontrada (x, y): {solution}")

🧰 Herramientas de Python: `NumPy` y `SciPy` (Uso Directo)
#

En la práctica, la función scipy.optimize.root es la herramienta estándar, ya que maneja internamente el ciclo de Newton.

from scipy.optimize import root

def fun_sistema_scipy(x):
    # x[0] = x, x[1] = y
    f1 = x[0]**2 + x[1]**2 - 4
    f2 = np.exp(x[0]) - x[1] - 1
    return np.array([f1, f2]) 

x_guess = [1.0, 1.0]
# El método 'hybr' (basado en Newton) calcula el Jacobiano numéricamente por defecto.
solucion = root(fun_sistema_scipy, x_guess, method='hybr') 

print(f"Raíz encontrada con root: {solucion.x}")

🌍 Ejemplo Práctico: Equilibrio de Reacciones
#

Problema
#

Encontrar el punto de equilibrio donde dos fuerzas se intersectan: $f_1(x, y) = x^3 - y^2 + 1 = 0$ y $f_2(x, y) = x y + x - 2 = 0$.

Solución con `SciPy.optimize.root` (Jacobiano Analítico)
#

Proporcionar el Jacobiano analítico es la práctica más robusta para optimizar velocidad y precisión.

Matriz Jacobiana Analítica:

$$\mathbf{J}(x, y) = \begin{pmatrix} 3x^2 & -2y \\ y+1 & x \end{pmatrix}$$

Código y Solución:

import numpy as np
from scipy.optimize import root

def fun_equilibrio(x):
    f1 = x[0]**3 - x[1]**2 + 1
    f2 = x[0] * x[1] + x[0] - 2
    return np.array([f1, f2])

def jac_equilibrio(x):
    J11 = 3 * x[0]**2
    J12 = -2 * x[1]
    J21 = x[1] + 1
    J22 = x[0]
    return np.array([[J11, J12], [J21, J22]])

x_guess = [1.0, 1.0] 

solucion = root(fun_equilibrio, x_guess, jac=jac_equilibrio, method='hybr')

print(f"--- Solución del Sistema No Lineal ---")
print(f"Raíz x: {solucion.x[0]:.6f}")
print(f"Raíz y: {solucion.x[1]:.6f}")
print(f"Residuos (F(x)): {solucion.fun}")
# Resultado: x ≈ 1.077209, y ≈ 1.340798

📝 Ejercicios Propuestos
#

A. Conceptuales y Fundamentos Teóricos
#

Definición del Jacobiano:
- Pregunta: Explica qué representa físicamente cada fila de la Matriz Jacobiana si el sistema $\mathbf{F}(\mathbf{x}) = \mathbf{0}$ proviene de un problema de equilibrio de fuerzas.
Ventaja Cuadrática:
- Pregunta: ¿Qué implica la convergencia cuadrática de Newton-Raphson? Si el error en la iteración $k$ es $0.01$, ¿cuál es el orden de magnitud esperado para el error en la iteración $k+1$?
Alternativa al Jacobiano Inverso:
- Pregunta: ¿Por qué las implementaciones modernas prefieren resolver un sistema lineal ($\mathbf{J} \boldsymbol{\Delta}\mathbf{x} = -\mathbf{F}$) en lugar de calcular explícitamente la inversa $\mathbf{J}^{-1}$?
Cálculo Analítico del Jacobiano:
- Tarea: Escribe la Matriz Jacobiana analítica $\mathbf{J}(x, y)$ para el sistema: $f_1(x, y) = \sin(x) + y^2 - 1 = 0$ y $f_2(x, y) = x^2 y - 0.5 = 0$.

B. Implementación con `SciPy` y `NumPy`
#

Sistema $3 \times 3$ Básico: Resuelve el siguiente sistema:
$$\begin{cases} f_1: 2x - y^2 + z = 0 \\ f_2: x^2 - y + 3z = 0 \\ f_3: x y z - 1 = 0 \end{cases}$$
- Tarea: Define la función residual y utiliza scipy.optimize.root con $\mathbf{x}^{(0)} = (1, 1, 1)$.
Punto de Intersección de Círculo y Parábola: Encuentra el punto de intersección (en el primer cuadrante) entre el círculo $x^2 + y^2 = 9$ y la parábola $y = x^2 - 1$.
- Tarea: Define las funciones residuales y usa root para encontrar la solución.
Newton Generalizado Asistido (Paso Único): Resuelve $f_1 = x^2 + y^2 - 4 = 0$ y $f_2 = x^3 - y - 1 = 0$ con $\mathbf{x}^{(0)}=(1, 1)$.
- Tarea:
  - Paso 1: Define $F(x, y)$ y $J(x, y)$ analíticamente.
  - Paso 2: Calcula $\mathbf{F}$ y $\mathbf{J}$ en $\mathbf{x}^{(0)}$.
  - Paso 3: Utiliza numpy.linalg.solve para obtener la corrección $\boldsymbol{\Delta}\mathbf{x}$.
  - Paso 4: Calcula $\mathbf{x}^{(1)} = \mathbf{x}^{(0)} + \boldsymbol{\Delta}\mathbf{x}$.
Importancia del Jacobiano Analítico:
- Tarea: Resuelve el ejercicio 7 de nuevo con scipy.optimize.root, pero esta vez proporciona el Jacobiano analítico (jac=jac_funcion).
- Pregunta: ¿Cómo se podría medir la diferencia en el tiempo de ejecución si se usara un sistema $100 \times 100$ (usando Jacobiano analítico vs. numérico)?
Análisis de Convergencia (Múltiples Raíces): El sistema $f_1(x, y) = x^2 + 4y^2 - 4 = 0$ y $f_2(x, y) = x - y + 1 = 0$ tiene dos soluciones reales.
- Tarea: Usa scipy.optimize.root para encontrar ambas soluciones, utilizando dos suposiciones iniciales muy diferentes: $\mathbf{x}^{(0)}=(0, 0)$ y $\mathbf{x}^{(0)}=(-2, 2)$.
Problemas Rígidos (Stiff) No Lineales:
- Pregunta de Análisis: ¿Por qué la resolución de los métodos Implícitos para EDOs (ej. Euler Implícito) a menudo requiere la aplicación del método de Newton-Raphson Generalizado en cada paso de tiempo? (Relaciona la EDO implícita con una ecuación algebraica no lineal para $y_{i+1}$).