Come ottenere una cancellazione più accurata
Cercherò di arrivare al punto, quindi fammi sapere se è rimasto qualcosa e hai bisogno di maggiori dettagli.
Sto risolvendo un paio di equazioni che non sono accoppiate esplicitamente , ma le loro corrispondenti variabili sconosciute, diciamo$x$ e $y$ deve soddisfare un'equazione differenziale:
$\dot x = x + y,$
dove i punti denotano derivata rispetto a una variabile indipendente, diciamo $t$.
L'equazione per $x$ è di secondo ordine, quindi si ottiene $x$ e $\dot x$da esso, e si può verificare se l'equazione di cui sopra è costantemente soddisfatta. Tuttavia (vedi la trama allegata), risulta che, indipendentemente dal metodo di integrazione che utilizzo da SciPy (quelli già implementati), l'uguaglianza di cui sopra smette di essere soddisfatta a un certo punto. Ciò è dovuto al fatto che$x$ e $y$ si cancellano a vicenda fino ad una precisione molto elevata, che sembra non essere raggiunta da nessuno dei metodi forniti da SciPy (l'ho verificato prendendo tutti i metodi e abbassando il più possibile la tolleranza assoluta e relativa. Nella trama allegata , il metodo impiegato è DOP853, che dovrebbe essere molto utile quando sono richieste tolleranze molto basse).
La mia domanda è se conosci un modo per migliorare l'accuratezza in modo che la cancellazione diventi più precisa (vorrei che l'equazione fosse convenientemente soddisfatta durante l'intero calcolo). Gli unici parametri che ho modificato finora sono state le tolleranze relative e assolute (e ovviamente i diversi metodi a nostra disposizione). C'è qualche parametro che mi manca e che potrebbe essere utile per questo?
Risposte
Non sono sicuro che ciò sia possibile con le librerie Python poiché stanno usando Fortran sotto il cofano e che non può essere facilmente ricompilato, ma la compilazione JIT Julia DifferentialEquations.jl specializza i solutori in base ai tipi di numero che gli dai. Ecco una dimostrazione di alcuni tipi strani come numeri razionali, MPFR BigFloats e ArbFloats (basati sulla libreria Arb) .
Puoi vederlo in azione nei grafici di convergenza di Feagin che dimostrano l'accuratezza di un metodo del 14 ° ordine$10^{-50}$tramite BigFloats . In Julia con BigFloats o ArbFloats, puoi setprecisionmodificare la precisione dei tipi di numero per ottenere la precisione di cui hai bisogno.
Mentre i metodi Julia sono molto veloci rispetto a SciPy (ordine e mezzo di grandezza) , e anche se si compilano per specializzarsi sui tipi di input, quindi viene generato uno speciale codice di ottimizzazione per il caso ad alta precisione, l'aritmetica ad alta precisione è ancora piuttosto costosa e dovresti tenerlo a mente. Specializzare l'integratore per questa gamma di alta precisione sarà piuttosto importante. Nota che se vuoi farlo, probabilmente consiglierei Vern9o uno dei metodi di estrapolazione multithread come ExtrapolationMidpointDeuflhard(eseguiranno il multithread tra le fchiamate che sarà più importante al diminuire della tolleranza), o forse il nuovo integratore simplettico del 16 ° ordineIRKGL16 .
Inoltre, se è necessaria un'aritmetica convalidata, è possibile utilizzare TaylorIntegration.jl per metodi Taylor di ordine elevato con limiti di precisione in virgola mobile sulla soluzione.