Cómo obtener una cancelación más precisa

Question

Intentaré ir al grano, así que avíseme si queda algo y necesita más detalles.

Estoy resolviendo un par de ecuaciones que no están acopladas explícitamente , pero sus correspondientes variables desconocidas, digamos$x$ y $y$ debe satisfacer una ecuación diferencial:

$\dot x = x + y,$

donde los puntos denotan derivada con respecto a una variable independiente, digamos $t$.

La ecuación para $x$ es de segundo orden, por lo que se obtiene $x$ y $\dot x$a partir de él, y se puede comprobar si la ecuación anterior se satisface de manera consistente. Sin embargo (ver el diagrama adjunto), resulta que, no importa qué método de integración use de SciPy (los que ya están implementados), la igualdad anterior deja de cumplirse en algún momento. Esto se debe al hecho de que$x$ y $y$ se cancelan entre sí con una precisión muy alta, que parece no lograrse con ninguno de los métodos proporcionados por SciPy (he comprobado esto tomando todos los métodos y reduciendo la tolerancia absoluta y relativa tanto como sea posible. En el gráfico adjunto , el método empleado es DOP853, que se supone que es muy útil cuando se requieren tolerancias muy bajas).

Mi pregunta es si conoce alguna forma de mejorar la precisión para que la cancelación sea más precisa (me gustaría que la ecuación se cumpliera convenientemente durante todo el cálculo). Los únicos parámetros que he cambiado hasta ahora fueron las tolerancias relativas y absolutas (y por supuesto los diferentes métodos a nuestra disposición). ¿Hay algún parámetro que me falte y que pueda ser útil para eso?

Answer 1

No estoy seguro de que esto sea posible con las bibliotecas de Python ya que están usando Fortran bajo el capó y eso no se puede recompilar fácilmente, pero la compilación JIT de Julia DifferentialEquations.jl especializa los solucionadores en función de los tipos de números que le da. Aquí hay una demostración de algunos tipos extraños como números racionales, MPFR BigFloats y ArbFloats (basados ​​en la biblioteca Arb) .

Puede ver esto en acción en los gráficos de convergencia de Feagin que demuestran la precisión de un método de orden 14 para$10^{-50}$a través de BigFloats . En Julia con BigFloats o ArbFloats, puede setprecisioncambiar la precisión de los tipos de números para obtener la precisión que necesita.

Si bien los métodos de Julia son muy rápidos en comparación con SciPy (orden y medio de magnitud) , y aunque se compilan para especializarse en tipos de entrada, por lo que se genera un código de optimización especial para el caso de alta precisión, la aritmética de alta precisión sigue siendo bastante costosa. y debes tener esto en cuenta. Especializar al integrador para este rango de alta precisión será bastante importante. Tenga en cuenta que si desea hacer esto, probablemente recomendaría Vern9uno de los métodos de extrapolación multiproceso como ExtrapolationMidpointDeuflhard(serán multiproceso entre fllamadas, lo que será más importante a medida que disminuya la tolerancia), o tal vez el nuevo integrador simpléctico de orden 16IRKGL16 .

Además, si necesita aritmética validada, puede usar TaylorIntegration.jl para métodos Taylor de alto orden con límites de precisión de punto flotante en la solución.