1. Introducción

Las ecuaciones de Maxwell describen el comportamiento de los campos eléctricos y magnéticos. En la práctica de simulación electromagnética, la mayoría de los métodos numéricos (FEM, FDTD, MME) convierten el problema continuo en un sistema lineal de la forma Ax = b. Este artículo muestra paso a paso cómo construir y resolver ese sistema usando herramientas de álgebra lineal en Python.

2. Fundamentos de las ecuaciones de Maxwell

• Gauss para el eléctrico: \(\nabla\cdot\mathbf{D}=\rho\)
• Gauss para el magnético: \(\nabla\cdot\mathbf{B}=0\)
• Faraday: \(\nabla\times\mathbf{E}= -\frac{\partial \mathbf{B}}{\partial t}\)
• Ampère‑Maxwell: \(\nabla\times\mathbf{H}= \mathbf{J}+\frac{\partial \mathbf{D}}{\partial t}\)

Al discretizar el dominio (malla de elementos o celda FDTD) cada ecuación se convierte en una relación algebraica entre nodos o celdas, generando matrices esparcidas extremadamente grandes.

4. Algoritmos de solución lineal

Dependiendo del tamaño y la densidad de A, se eligen diferentes métodos:

• Directos (LU, Cholesky): Ideales para matrices < 10⁶ elementos, garantizan solución exacta (hasta error numérico).
• Iterativos (Conjugate Gradient, GMRES): Recomendados para matrices esparcidas >10⁶, requieren pre‑condicionadores.
• Eigenvalue solvers: Cuando se busca modos propios (cavidades resonantes).

Python ofrece scipy.linalg (directo) y scipy.sparse.linalg (iterativo).

5. Implementación en Python

A continuación se muestra una plantilla genérica que puede adaptarse a FEM o FDTD.

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
# --------------------------------------------------
# 1. Construcción de la matriz esparcida K
# --------------------------------------------------
# Ejemplo simple 2‑D usando elementos bilineales en una cuadrícula
N = 100  # número de nodos por dirección
h = 1.0 / (N-1)
size = N * N
# Operador laplaciano 5‑point stencil
main_diag = -4 * np.ones(size)
side_diag = np.ones(size-1)
up_down_diag = np.ones(size-N)
# Corregir los cruces de fila
for i in range(1, N):
    side_diag[i*N-1] = 0
K = sp.diags([main_diag, side_diag, side_diag, up_down_diag, up_down_diag],
             [0, -1, 1, -N, N], format='csr') / h**2
# --------------------------------------------------
# 2. Fuente (vector b)
# --------------------------------------------------
b = np.zeros(size)
# Fuente puntual en el centro
i_center = N//2 * N + N//2
b[i_center] = 1.0
# --------------------------------------------------
# 3. Solución directa (LU) – para pruebas pequeñas
# --------------------------------------------------
phi_direct = spla.spsolve(K, b)
# --------------------------------------------------
# 4. Solución iterativa (CG) con pre‑condicionador Incomplete LU
# --------------------------------------------------
M = spla.spilu(K)               # factorización ILU
M_x = lambda x: M.solve(x)     # función de aplicación
M_inv = spla.LinearOperator((size, size), matvec=M_x)
phi_cg, info = spla.cg(K, b, M=M_inv, tol=1e-8, maxiter=1000)
print('CG convergió en', info, 'iteraciones')
# --------------------------------------------------
# 5. Visualización rápida (requiere matplotlib)
# --------------------------------------------------
import matplotlib.pyplot as plt
phi_grid = phi_cg.reshape((N, N))
plt.imshow(phi_grid, origin='lower', cmap='viridis')
plt.colorbar(label='Potencial')
plt.title('Solución de Poisson (analógica a Maxwell)')
plt.show()

El fragmento anterior resuelve el problema de Poisson, que es la forma escalar de la ecuación de Helmholtz que aparece al eliminar el tiempo en Maxwell bajo condiciones armónicas. Cambiando K por la matriz curl‑curl y añadiendo la constante de onda k, se obtiene la solución de frecuencia de los campos.

7. Buenas prácticas y optimización

• Usar tipos esparcidos desde el inicio: evita conversiones costosas de numpy.ndarray a scipy.sparse.
• Escalar la matriz (pre‑escalado) para mejorar la condición numérica.
• Elección del pre‑condicionador: ILU funciona bien para problemas de Poisson; para Helmholtz se recomiendan AMG o HYPRE.
• Vectorización y uso de BLAS: las operaciones de producto matriz‑vector son críticas; verifica que NumPy esté vinculada a MKL/OpenBLAS.
• Persistencia de la factorización: si resuelves varias cargas con la misma K, guarda la factorización LU/ILU y reutilízala.

8. Escalabilidad y paralelismo

Para simulaciones de gran escala (>10⁷ DOF) se recomienda:

• Distribuir la matriz con mpi4py y usar solvers iterativos paralelos como petsc4py o pyamg con MPI.
• Ejecutar en GPU mediante cupy o torch.linalg para productos Ax masivos.
• Utilizar algoritmos de dominio‑descompuesto (Schwarz, FETI) cuando la malla es heterogénea.

Ejemplo rápido de mpi4py con CG:

from mpi4py import MPI
import petsc4py
petsc4py.init()
from petsc4py import PETSc
comm = MPI.COMM_WORLD
A = PETSc.Mat().createAIJ([N,N], comm=comm)  # ensamblado distribuido
b = PETSc.Vec().create(comm=comm)
# ... ensamblar A y b ...
ksp = PETSc.KSP().create(comm=comm)
ksp.setOperators(A)
ksp.setType('cg')
ksp.getPC().setType('hypre')
solution = b.duplicate()
ksp.solve(b, solution)