Pesos que salvan el gradiente en mallas no estructuradas — aproximación por mínimos cuadrados móviles (MLS)

El día en que un gradiente se rompió en una malla no estructurada#

A las dos de la tarde, sobre una malla tetraédrica no estructurada, el gradiente LSQ (mínimos cuadrados) saltó negativo en una sola celda. La vecina marcaba 0.3, la siguiente 0.5, pero la del medio caía a -4.1. Al revisar la malla, los siete vecinos de esa celda estaban a distancias muy desiguales. Un punto cercano entraba en las ecuaciones normales con el mismo peso que seis lejanos.

La solución ya existía desde 1981, con otro nombre y en otra disciplina.

Al mover h, el orden del polinomio o el número de muestras N, el RMSE de la derivada se mueve con ellos.

Reducir $h$ a 0.025 deja a algunos puntos de evaluación con demasiado pocos nodos activos: $\mathbf{A}$ se vuelve casi singular y el RMSE explota. Subir $h$ a 0.20 promedia la amplitud de la función y embota la derivada. Entre ambos extremos hay un sweet spot estrecho.

1981 — un préstamo del ajuste de superficies#

Lancaster y Salkauskas hicieron una promesa en un artículo sobre reconstrucción de superficies: "en cada punto de evaluación, reajustar un polinomio usando solo los datos cercanos a ese punto". Lo llamaron moving least-squares, MLS. Era una herramienta de gráficos para construir superficies suaves.

Veinte años después el EFG (Element-Free Galerkin) de Belytschko tomó la misma ecuación como esqueleto del análisis estructural sin malla. Diez años más y el FVM no estructurado se la llevó para reconstruir gradientes. La derivada con corrección de consistencia de SPH vive en el mismo lugar. El punto de partida común — puntos en lugar de mallas — invoca la misma ecuación en cada campo.

Una función de peso transforma distancia en influencia#

La ecuación de partida del MLS es simple. Cerca del punto de evaluación $\mathbf{x}^*$, se ajusta la función con una base polinómica.

$\mathbf{p}$ es el vector de la base polinómica (lineal en 2D: $[1, \xi, \eta]$; cuadrática en 2D: $[1, \xi, \eta, \xi^2/2, \xi\eta, \eta^2/2]$), y $\mathbf{a}$ es el vector de coeficientes.

Lo crucial: $\mathbf{a}$ es una función del punto de evaluación $\mathbf{x}^*$. Si el punto se mueve, los coeficientes se vuelven a resolver. De ahí el "moving" en mínimos cuadrados móviles.

Los coeficientes se obtienen minimizando la suma ponderada de residuos al cuadrado, $J$.

La función de peso $w$ hace el trabajo pesado. Los puntos cercanos a $\mathbf{x}^*$ reciben pesos grandes; los lejanos, casi cero. Gaussiana $w = \exp(-(r/h)^2)$, spline cúbica, polinomio de cuarto orden — todas comparten un rasgo. Más allá de una distancia $h$ se cortan a cero. Esa $h$ es el kernel radius (radio de soporte).

Ecuaciones normales y la aparición de la función de forma#

La condición $\partial J / \partial \mathbf{a} = 0$ produce las ecuaciones normales:

con

$\mathbf{u}$ contiene los valores nodales y $w_i = w(\|\mathbf{x}^* - \mathbf{x}_i\|)$. Sustituyendo los coeficientes:

$\Phi_i(\mathbf{x}^*)$ es la función de forma (shape function). Trae una propiedad pequeña pero molesta: $\Phi_i(\mathbf{x}_j) \neq \delta_{ij}$. La interpolación no es exacta en los nodos. Eso complica al EFG cuando se imponen condiciones Dirichlet. Para la reconstrucción de gradientes en CFD no es relevante: lo que queremos es una derivada suave, no una interpolación nodal exacta.

Python — derivar el seno desde 25 puntos#

Eligiendo la base en forma de Taylor, $p_k(\xi) = \xi^k/k!$, los coeficientes pasan a ser las derivadas directamente — $a_k \approx u^{(k)}(\mathbf{x}^*)$.

import numpy as np
from math import factorial
 
def mls_taylor_1d(xs, us, x_star, h, order=2):
    """Mínimos cuadrados ponderados con base de Taylor en x_star.
    El a[k] devuelto aproxima u^(k)(x_star).
    """
    m = order + 1
    A = np.zeros((m, m))
    B = np.zeros(m)
    for x_i, u_i in zip(xs, us):
        dx = x_i - x_star
        w = np.exp(-(dx / h) ** 2)
        if w < 1e-10:
            continue
        p = np.array([(dx ** k) / factorial(k) for k in range(m)])
        A += w * np.outer(p, p)
        B += w * p * u_i
    return np.linalg.solve(A, B)
 
# Validación — recuperar la derivada de sin(2 pi x)
np.random.seed(7)
xs = np.sort(np.random.uniform(0, 1, 25))
us = np.sin(2 * np.pi * xs)
 
x_eval = np.linspace(0.05, 0.95, 200)
du_hat = np.array([mls_taylor_1d(xs, us, xq, h=0.07, order=2)[1] for xq in x_eval])
du_true = 2 * np.pi * np.cos(2 * np.pi * x_eval)
 
rmse = np.sqrt(np.mean((du_hat - du_true) ** 2))
print(f"derivative RMSE: {rmse:.3e}")
# derivative RMSE: 2.1e-01   (h=0.07, order=2)

Veinticinco puntos aleatorios, RMSE de la derivada alrededor de 0.21. La amplitud es $2\pi \approx 6.28$, así que el error relativo ronda el 3%. Subiendo order a 3 y $h$ a 0.10 baja todavía más. Pero si $h$ se agranda demasiado, aparece otro problema.

En la gaussiana $w=\exp(-(r/h)^{2})$ el peso cae a la mitad donde $(r/h)^{2}=\ln 2$. Con $h=0.07$ ese radio vale unos 0.058. El espaciado medio de los 25 puntos ronda 0.042. El fragmento descarta un punto cuando $w<10^{-10}$, a unos $4.8h$: 0.34 si $h=0.07$, y 0.12 si $h=0.025$.

Orden 2 significa tres incógnitas. La ventana de $h=0.025$ mide unos 0.24 de lado a lado, cinco o seis puntos con un espaciado de 0.042. La semilla 7 deja huecos bajo ese promedio, y ahí $\mathbf{A}$ se pone mal condicionada. Dividir el RMSE 0.21 de la derivada entre la amplitud $2\pi\approx 6.28$ da 0.033, cerca del 3%.

Tomar $h$ como 1.5 veces el espaciado medio da $0.042\times 1.5\approx 0.063$, casi el 0.07 del fragmento. Llevar $h$ a 3 veces, unos 0.13, hace que el radio de medio peso pase dos vecinos y embote la derivada alta.

El caso sin peso figura en Green–Gauss y mínimos cuadrados ponderados. El −4.1 de esa tarde fue el mismo coeficiente en un punto cercano y uno lejano.

El radio de núcleo h — qué destruye cada extremo#

En reconstrucción de gradientes, $h$ negocia dos cosas.

• $h$ pequeño: el número de condición de $\mathbf{A}$ se dispara. Uno o dos puntos cercanos cargan todo el peso, y se intenta ajustar un polinomio de grado $m$ con menos de $m$ puntos activos. El resultado se vuelve sensible al ruido.
• $h$ grande: la información de puntos lejanos se mezcla y promedia. El momento de orden cero sobrevive, pero las derivadas de orden superior se embotan.

En la práctica, se elige $h$ entre 1.5 y 3 veces el espaciado nodal local medio, de modo que el soporte contenga entre 1.5 y 2 veces $m$ nodos activos. En mallas no estructuradas hay que adaptar $h$ celda a celda — un $h$ global casi siempre falla.

Vamos al caso 2D. Veintiocho nodos dispersos llevan $u = \sin(2\pi x)\cos(2\pi y)$. Arrastra el punto amarillo de evaluación y observa cómo cambian los pesos.

Si el punto se mueve a la esquina superior izquierda vacía, el número de nodos activos cae por debajo de $m$ y $\mathbf{A}$ se vuelve singular. Llevarlo a una zona densa devuelve el error a valores pequeños. Esta es la debilidad estructural del MLS en mallas no estructuradas, y la razón por la que en la práctica un $h$ adaptativo no es opcional.

Ecos en SPH, EFG y FVM#

El MLS vive hoy en tres lugares bajo nombres distintos.

• Corrección de consistencia en SPH: la derivada estándar de SPH pierde incluso la precisión de primer orden con distribuciones no uniformes de partículas. SPH corregido y reproducing-kernel SPH inyectan una corrección por mínimos cuadrados ponderados para restaurar la consistencia de primer orden — esencialmente MLS de orden uno.
• EFG / MLPG, análisis estructural sin malla: la función de forma es literalmente MLS. La no exactitud en la interpolación nodal complica las condiciones Dirichlet; multiplicadores de Lagrange o métodos de penalización son la solución habitual.
• Reconstrucción de gradientes en FVM no estructurado: Compact Least-Squares Reconstruction resuelve las ecuaciones normales sobre los centros de celda y extrae no solo gradientes sino también curvatura. La opción leastSquares de OpenFOAM es una variante de MLS con base polinómica de orden uno.

Tres campos, una sola familia de ecuaciones normales. Solo cambian la función de peso y la dimensión de la base polinómica.

Volviendo al gradiente roto#

De vuelta a la celda inicial. El gradiente LSQ saltó a -4.1 porque dos vecinos cercanos entraban en las ecuaciones normales con el mismo peso que cinco lejanos. Con un peso por distancia, -4.1 pasó a 0.4 y la divergencia desapareció.

Localidad por función de peso, ecuaciones normales reformuladas en cada punto y la no separabilidad de la función de forma — tres piezas que el MLS prestó desde los gráficos. Hoy toda la integración sin malla y la mitad de la integración no estructurada se apoyan en ese andamiaje.

La celda cuyo gradiente salta a −4.1 cuando cada vecino comparte un peso#

Dos vecinos cercanos y cinco lejanos entraron con el mismo coeficiente, y el gradiente de esa celda saltó a −4.1. Un peso por distancia lo devolvió a 0.4. Si $h$ queda muy por debajo del espaciado nodal, $\mathbf{A}$ se inclina a singular. El signo puede verse bien mientras la magnitud vuelve a temblar.

Con el gradiente de vuelta, la decisión siguiente es cuánto de esa pendiente puede entrar en la celda. El tramo en el que una reconstrucción de segundo orden fabrica extremos falsos en una malla no estructurada está en dónde corta un limitador de Barth–Jespersen. Devolver la pendiente con pesos y omitir el limitador hace que esa celda vuelva a crecer un extremo.