Numpy: Cálculo

NumPy proporciona un conjunto de herramientas poderosas para realizar cálculo numérico en Python. A diferencia del cálculo simbólico (realizado por bibliotecas como SymPy), el cálculo numérico se centra en aproximaciones computacionales de problemas matemáticos.

Las aplicaciones del cálculo numérico son extensas:

• Simulaciones físicas
• Análisis financiero
• Procesamiento de señales
• Modelado climático
• Mecánica de fluidos

El módulo numpy no incluye funciones específicas para cálculo diferencial e integral como parte de su API principal, pero proporciona las estructuras de datos y operaciones fundamentales para implementar estos métodos.

Diferenciación numérica: cálculo de derivadas

La derivada representa la tasa de cambio de una función. En NumPy, podemos aproximarla mediante diferencias finitas.

El método más sencillo es la diferencia hacia adelante:

import numpy as np

def derivada_adelante(f, x, h=1e-5):
    return (f(x + h) - f(x)) / h

# Ejemplo: derivada de sin(x) en x=π/4
x = np.pi/4
def f(x): return np.sin(x)
print(f"Derivada numérica: {derivada_adelante(f, x)}")
print(f"Valor exacto (cos(π/4)): {np.cos(x)}")

Una aproximación más precisa es la diferencia central:

def derivada_central(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

Para calcular derivadas de orden superior, podemos aplicar recursivamente el operador de diferenciación:

def segunda_derivada(f, x, h=1e-5):
    return (f(x + h) - 2 * f(x) + f(x - h)) / (h**2)

También podemos calcular derivadas parciales de funciones multivariables:

def derivada_parcial(f, x, y, h=1e-5, respecto_a='x'):
    if respecto_a == 'x':
        return (f(x + h, y) - f(x, y)) / h
    else:  # respecto a 'y'
        return (f(x, y + h) - f(x, y)) / h

Gradientes y derivadas de funciones vectoriales

El gradiente es una generalización de la derivada para funciones multivariables. Representa la dirección de máximo crecimiento de una función.

Podemos implementar una función para calcular gradientes aproximados:

def gradiente(f, punto, h=1e-5):
    """Calcula el gradiente de f en el punto dado."""
    n = len(punto)
    grad = np.zeros(n)
    
    for i in range(n):
        punto_adelante = punto.copy()
        punto_atras = punto.copy()
        punto_adelante[i] += h
        punto_atras[i] -= h
        grad[i] = (f(punto_adelante) - f(punto_atras)) / (2 * h)
    
    return grad

# Ejemplo: gradiente de f(x,y) = x²+y² en (1,2)
def f(punto):
    x, y = punto
    return x**2 + y**2

punto = np.array([1.0, 2.0])
print(f"Gradiente numérico: {gradiente(f, punto)}")
print(f"Gradiente exacto: [2x, 2y] = {[2*punto[0], 2*punto[1]]}")

Para funciones de una variable definidas sobre un array, podemos usar np.gradient:

# Calcular la derivada de sin(x) sobre un rango
x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x)
dy_dx = np.gradient(y, x)

# La derivada de sin(x) es cos(x)
import matplotlib.pyplot as plt
plt.plot(x, dy_dx, 'r-', label='Derivada numérica')
plt.plot(x, np.cos(x), 'b--', label='cos(x) exacto')
plt.legend()
plt.title('Comparación de derivada numérica vs exacta')

Integración numérica con NumPy

La integración numérica permite calcular aproximaciones de integrales definidas. El método más simple es la regla del trapecio.

def regla_trapecio(f, a, b, n=100):
    """Integra f de a a b usando n trapecios."""
    x = np.linspace(a, b, n+1)
    y = f(x)
    return (b - a) * (np.sum(y) - 0.5*y[0] - 0.5*y[-1]) / n

# Integrando sin(x) de 0 a pi
def f(x): return np.sin(x)
integral = regla_trapecio(f, 0, np.pi)
print(f"∫sin(x)dx de 0 a π ≈ {integral}")
print(f"Valor exacto: {2.0}")

NumPy proporciona la función np.trapz para aplicar esta regla de forma más directa:

x = np.linspace(0, np.pi, 1000)
y = np.sin(x)
integral = np.trapz(y, x)
print(f"Usando np.trapz: {integral}")

Para integraciones más precisas, podemos implementar la regla de Simpson:

def regla_simpson(f, a, b, n=100):
    """Integra f de a a b usando n intervalos con la regla de Simpson."""
    if n % 2 != 0:
        n += 1  # Asegurarse de que n sea par
    
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    y = f(x)
    
    resultado = y[0] + y[-1]
    resultado += 4 * np.sum(y[1:-1:2])  # índices impares
    resultado += 2 * np.sum(y[2:-1:2])  # índices pares
    
    return resultado * h / 3

# Ejemplo
integral_simpson = regla_simpson(f, 0, np.pi)
print(f"Regla de Simpson: {integral_simpson}")

Integrales dobles y múltiples

Las integrales múltiples pueden calcularse mediante la aplicación recursiva de la integración numérica:

def integral_doble(f, ax, bx, ay, by, nx=20, ny=20):
    """Calcula ∫∫f(x,y)dxdy sobre [ax,bx]×[ay,by]."""
    x = np.linspace(ax, bx, nx)
    y = np.linspace(ay, by, ny)
    hx = (bx - ax) / (nx - 1)
    hy = (by - ay) / (ny - 1)
    
    resultado = 0
    for i in range(nx):
        for j in range(ny):
            # Factor del peso del punto (x[i], y[j])
            peso = 1.0
            if i == 0 or i == nx-1: peso *= 0.5
            if j == 0 or j == ny-1: peso *= 0.5
            
            resultado += peso * f(x[i], y[j])
    
    return resultado * hx * hy

# Ejemplo: ∫∫(x+y)dxdy sobre [0,1]×[0,1]
def f(x, y): return x + y
integral_2d = integral_doble(f, 0, 1, 0, 1)
print(f"Integral doble: {integral_2d}")
print(f"Valor exacto: 1.0")

Ecuaciones diferenciales con NumPy

Las ecuaciones diferenciales modelan sistemas dinámicos y procesos físicos. NumPy nos permite implementar métodos numéricos para resolverlas.

El método más sencillo es el método de Euler:

def euler_method(f, y0, t_range, h):
    """
    Resuelve una EDO dy/dt = f(t, y) con el método de Euler.
    f: función que define la derivada
    y0: condición inicial
    t_range: (t_inicio, t_fin)
    h: tamaño del paso
    """
    t_inicio, t_fin = t_range
    t = np.arange(t_inicio, t_fin, h)
    y = np.zeros_like(t)
    y[0] = y0
    
    for i in range(1, len(t)):
        y[i] = y[i-1] + h * f(t[i-1], y[i-1])
    
    return t, y

# Ejemplo: resolver dy/dt = -2y con y(0) = 1
# Solución exacta: y = e^(-2t)
def f(t, y): return -2 * y
t, y = euler_method(f, 1.0, (0, 2), 0.01)

# Comparar con solución exacta
y_exacta = np.exp(-2 * t)
error = np.abs(y - y_exacta)
print(f"Error máximo: {np.max(error)}")

Un método más preciso es Runge-Kutta de cuarto orden (RK4):

def runge_kutta4(f, y0, t_range, h):
    """Resuelve una EDO con el método RK4."""
    t_inicio, t_fin = t_range
    t = np.arange(t_inicio, t_fin, h)
    y = np.zeros_like(t)
    y[0] = y0
    
    for i in range(1, len(t)):
        ti = t[i-1]
        yi = y[i-1]
        
        k1 = f(ti, yi)
        k2 = f(ti + h/2, yi + h*k1/2)
        k3 = f(ti + h/2, yi + h*k2/2)
        k4 = f(ti + h, yi + h*k3)
        
        y[i] = yi + h * (k1 + 2*k2 + 2*k3 + k4) / 6
    
    return t, y

# Resolver la misma ecuación con RK4
t_rk4, y_rk4 = runge_kutta4(f, 1.0, (0, 2), 0.01)
error_rk4 = np.abs(y_rk4 - np.exp(-2 * t_rk4))
print(f"Error máximo RK4: {np.max(error_rk4)}")

Para sistemas de ecuaciones diferenciales, extendemos estos métodos trabajando con vectores:

def sistema_rk4(f, y0, t_range, h):
    """
    Resuelve un sistema de EDOs con RK4.
    f: función que retorna un vector con las derivadas
    y0: vector de condiciones iniciales
    """
    t_inicio, t_fin = t_range
    t = np.arange(t_inicio, t_fin, h)
    n = len(t)
    y = np.zeros((n, len(y0)))
    y[0] = y0
    
    for i in range(1, n):
        ti = t[i-1]
        yi = y[i-1]
        
        k1 = f(ti, yi)
        k2 = f(ti + h/2, yi + h*k1/2)
        k3 = f(ti + h/2, yi + h*k2/2)
        k4 = f(ti + h, yi + h*k3)
        
        y[i] = yi + h * (k1 + 2*k2 + 2*k3 + k4) / 6
    
    return t, y

# Ejemplo: sistema de un oscilador armónico
# dy1/dt = y2
# dy2/dt = -y1
def oscilador(t, y):
    dy = np.zeros(2)
    dy[0] = y[1]
    dy[1] = -y[0]
    return dy

# Condiciones iniciales: y1(0)=0, y2(0)=1
t, y = sistema_rk4(oscilador, np.array([0, 1]), (0, 10), 0.01)

# La solución exacta es y1 = sin(t), y2 = cos(t)

Aproximación de funciones y series

Las series de Taylor permiten aproximar funciones mediante polinomios. Veamos cómo implementarlas con NumPy:

def serie_taylor(f, derivadas, x0, x, n_terminos=5):
    """
    Aproxima f(x) usando n_terminos de la serie de Taylor en torno a x0.
    f: función a aproximar
    derivadas: lista de funciones con las derivadas de f
    """
    aproximacion = f(x0)
    factorial = 1
    
    for i in range(1, n_terminos):
        factorial *= i
        termino = derivadas[i-1](x0) * (x - x0)**i / factorial
        aproximacion += termino
    
    return aproximacion

# Ejemplo: aproximar sin(x) en torno a x=0
def f(x): return np.sin(x)
# Primeras derivadas de sin(x): cos(x), -sin(x), -cos(x), sin(x)
derivadas = [
    lambda x: np.cos(x),
    lambda x: -np.sin(x),
    lambda x: -np.cos(x),
    lambda x: np.sin(x)
]

x = np.linspace(-np.pi/2, np.pi/2, 100)
y_exacto = np.sin(x)
y_aprox = serie_taylor(f, derivadas, 0, x, 5)

print(f"Error máximo: {np.max(np.abs(y_exacto - y_aprox))}")

Otra forma de aproximar funciones es mediante interpolación polinómica:

def interpolacion_lagrange(x, y, x_eval):
    """
    Realiza interpolación polinómica de Lagrange.
    x, y: puntos de interpolación
    x_eval: puntos donde evaluar el polinomio
    """
    n = len(x)
    result = np.zeros_like(x_eval, dtype=float)
    
    for i in range(n):
        # Término L_i(x)
        term = y[i]
        for j in range(n):
            if i != j:
                term = term * (x_eval - x[j]) / (x[i] - x[j])
        result += term
    
    return result

# Ejemplo: interpolar sin(x) con 5 puntos
x_datos = np.linspace(0, np.pi, 5)
y_datos = np.sin(x_datos)

x_eva