Por qué las recurrencias de Bessel fallan en tu código

El problema fundamental: dos soluciones, una sola confiable

Una ecuación en diferencias lineal de segundo orden tiene dos soluciones independientes, igual que una ecuación diferencial de segundo orden. Cualquier solución particular es una combinación lineal de esas dos.

Aquí está la trampa: si una solución crece con n y la otra decae, lo que normalmente querés calcular es la que decae, la llamada solución mínima. Pero al hacerlo, el redondeo de punto flotante introduce una pequeña contaminación de la solución que crece. En cada paso, esa contaminación se amplifica. Tras suficientes iteraciones, domina el resultado y diverge del valor teórico.

Este fenómeno no es exclusivo de las funciones de Bessel. Aparece en cualquier recurrencia donde el cociente entre las dos soluciones independientes tiende a cero o a infinito en función de n. La estabilidad, por tanto, depende de la dirección en que se aplica la recurrencia.

El caso clásico: funciones de Bessel Jn y Yn

Cuando mirás un gráfico de funciones de Bessel en un libro de texto, ves que Jn(x) y Yn(x) se comportan aproximadamente igual en función de x, como seno y coseno. Y eso es cierto.

Pero no es cierto en función de n para un x fijo. A medida que n crece, Jn(x) decae a cero mientras que Yn(x) diverge hacia −∞.

Esa diferencia de comportamiento como función de n es exactamente la fuente de la inestabilidad numérica. La recurrencia hacia adelante (calcular valores de orden superior a partir de valores de orden inferior) es estable para Yn pero inestable para Jn. En dirección inversa, la situación se invierte: estable para Jn, inestable para Yn.

Por qué la dirección de la recurrencia importa

En una nota anterior del mismo autor se mostró empíricamente que la estabilidad depende tanto del tipo de función de Bessel como de la dirección en que se aplica la recurrencia. Este artículo cierra el círculo y explica el motivo.

• Dirección hacia adelante (de orden bajo a orden alto): amplifica cualquier componente de la solución que crece con n. Si la solución que querés es la que decae (como J_n), tu resultado se contamina rápidamente y termina dominado por ruido.
• Dirección inversa (de orden alto a orden bajo): amplifica cualquier componente de la solución que crece hacia atrás. Si la solución que querés es la que crece hacia adelante (como Y_n), terminás amplificando el componente equivocado.

La regla práctica: aplicá la recurrencia hacia la solución mínima, no en dirección que se aleja de ella.

Miller's algorithm y otras técnicas

Existen técnicas numéricas para calcular la solución mínima (la que no diverge) de forma estable. Una mencionada en el artículo es el algoritmo de Miller: aplica la recurrencia en la dirección que sí es estable y luego renormaliza el resultado para imponer la condición de contorno correcta. La idea es generalizable a cualquier recurrencia con el mismo patrón de estabilidad direccional.

El principio general que deja el artículo: si el cociente de las dos soluciones independientes tiende a cero o infinito en función de n, esperá problemas de estabilidad. La única forma robusta de detectarlos es comparar empíricamente resultados en distintas direcciones y precisiones antes de confiar en el cálculo.

Qué significa esto para tu startup

Aunque pocas startups calculan funciones de Bessel a mano, el principio detrás del artículo aparece en tres lugares comunes del software de negocio:

• Cálculos financieros y de riesgo: sumar muchos términos pequeños con magnitudes muy distintas amplifica errores de redondeo. Si tu SaaS hace pricing dinámico, simulación de Monte Carlo o cálculo de provisiones, vale la pena saber cómo se está propagando el error en cada paso.
• Modelos de ML iterativos: algunos algoritmos (factorización de matrices, optimizadores con momentum, simulaciones físicas) iteran muchas veces sobre los mismos valores. Errores pequeños se acumulan de forma matemáticamente análoga a una recurrencia mal condicionada.
• Integración con librerías científicas: si usás NumPy, SciPy o similares, las funciones especiales ya están implementadas con el algoritmo estable por defecto. Pero las versiones optimizadas para velocidad a veces intercambian precisión por performance, así que vale verificar en el rango de valores que tu código realmente procesa.

Dos acciones concretas que podés tomar hoy:

1. Auditá tus pipelines numéricos buscando operaciones donde sumás muchos términos de magnitudes muy distintas. Si encontrás, considerá usar doble precisión (float64) en vez de float32, o técnicas como Kahan summation para reducir el error acumulado sin cambiar de tipo de dato.
2. Si dependés de una librería para funciones especiales (Bessel, gamma, erf, etc.), verificá que esté usando el algoritmo estable para el rango de valores que tu código procesa. Las librerías maduras lo hacen por defecto; las optimizadas para velocidad, no siempre.

Fuentes

• Numerical (In)Stability of Recurrence Relations