Ok, I might try that with interval, to see what I can get. I’m not sure of the math of Richardson extrapolation. I’m thinking it must use standard analysis to find the limit. I was using an infentesimal to find a one sided limit.
using SymEngine
ε=10^-308 #smallest 64 bit number
subs(x/((x^2)-1),x=>ε)
then I get Inf.0 as a result (technically incorrect, although if I could get it to say ω.0 or H.0, then it be a correct hyperreal answer.
I suppose I could do the same with extrapolate.
Richardson extrapolation works fine for this function, giving the correct result (zero):
julia> extrapolate(x -> x / (x^2 - 1), 0.1, x0=0)
(0.0, 0.0)
These course notes are a pretty good explanation of how it works.
I’m not sure why you are using 1e-308 (which is between floatmin(Float64) ≈ 2.22507e-308 and nextfloat(0.0) = 5.0e-324). It is probably not what you want — in non-symbolic calculations it might lead to erroneous results due to underflow or similar roundoff errors. For example, consider (e^x - 1) / x, which goes to 1 as x \to 0. If you plug in 1e-308 (or any sufficiently small number) you get the wrong answer (0.0) due to rounding:
julia> f(x) = (exp(x) - 1) / x
f (generic function with 1 method)
julia> f(1e-308)
0.0
julia> f(1e-20)
0.0
whereas Richardson extrapolation gives the correct answer:
julia> extrapolate(f, 0.1, x0=0)
(1.0000000000021632, 5.938805003324887e-12)
(In this particular case, you could also use the expm1 function to compute e^x - 1 more accurately for small x, but of course such a special function is not always available.)
Ok thanks, I’m not familiar with expm1, but I can do this. Is underflow because the number is too large, I was afraid there would be problems with a 256 bit infinitesimal, but I will definitely look into this extrapolation function, it seems like it would work better.
I was going to say the limit is infinity, but the right limit of x=1 is infinity, I think I can put in Inf, and get 0, but if I get the domain, then this would be more useful.