Hello, I’m a new user of Julia.

I want to invert numerically a function and to compute the inverse on a ‘large’ array.

For the moment I studied only two options: using a Newton-Rhapson algorithm (I know the first derivative analytically) and using interpolation.

The function I have to invert is monotonically increasing and it is called `luminosity_distance_at_z`

in the example below.

```
using BenchmarkTools
using Interpolations
using QuadGK
function NewtonRhapson(f, df, y0; x_guess=1.0, tol=1e-6, maxIter = 1000)
# found x0 : f(x0) = y0
ff(x) = f(x) - y0
x = x_guess
ffx = ff(x)
iter = 0
while abs(ffx) > tol && iter < maxIter
x -= ffx/df(x) # Iteration
ffx = ff(x) # Precompute ff(x)
iter += 1
end
return x
end
function inverse_interpolation(f, y0s)
xs = vcat( range(start = 0, stop = 9.99e-9, length = 125),
10 .^(range(start = -8, stop = log10(7.99), length = 1000)),
10 .^(range(start = log10(8), stop = 8, length = 125)) )
ys = f.(xs)
inv = interpolate((ys,), xs, Gridded(Linear()))
return inv.(y0s)
end
function E_at_z(z::Real)
sqrt(0.3*(1+z)^3 + 0.69991 + 9e-5*(1+z)^4)
end
function integral_of_E_inv(z::Real)
QuadGK.quadgk(z -> 1/E_at_z(z), 0., z)[1]
end
const c_light = 299792.458 # km/s
function comoving_radial_distance_at_z(z::Real)
dH = c_light/67
return dH * integral_of_E_inv(z)
end
function differential_luminosity_distance_at_z(z::Real)
dH = c_light/67
comoving_radial_distance_at_z(z) + (dH/E_at_z(z))*(1+z)
end
function luminosity_distance_at_z(z::Real)
comoving_radial_distance_at_z(z)*(1+z)
end
```

When I have to evaluate the inverse function on a single point, the NR algorithm wins.

When I have to evaluate the inverse on a ‘large’ set of points, the interpolation algorithm seems to perform better

I think that the difference is due to the fact that the interpolation method computes `f(x)`

1250 times (in this example), while the NR algorithm computes `f(x)`

two times each step needed to evaluate each point of the array (and a few steps per points are usually needed).

Since the NR is more precise, are there any ways to make it as fast as the interpolation one when evaluated on a large set of points?

Since I’m new I thought that maybe I’m using wrongly the vectorization of the NR algorithm via `broadcast`

and there could be more efficient ways to vectorize it.