Is there any way to bound the region searched by NLsolve in Julia?

Question

I'm trying to find one of the roots of a nonlinear (roughly quartic) equation. The equation always has four roots, a pair of them close to zero, a large positive, and a large negative root. I'd like to identify either of the near zero roots, but nlsolve, even with an initial guess very close to these roots, seems to always converge on the large positive or negative root.

A plot of the function essentially looks like a constant negative value, with a (very narrow) even-ordered pole near zero, and gradually rising to cross zero at the large positive and negative roots.

Is there any way I can limit the region searched by nlsolve, or do something to make it more sensitive to the presence of this pole in my function?

EDIT: Here's some example code reproducing the problem:

using NLsolve

function f!(F,x)
    x = x[1]
    F[1] = -15000 + x^4 / (x+1e-5)^2
end
# nlsolve will find the root at -122
nlsolve(f!,[0.0])

As output, I get:

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.0]
 * Zero: [-122.47447713915808]
 * Inf-norm of residuals: 0.000000
 * Iterations: 15
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 16
 * Jacobian Calls (df/dx): 6

We can find the exact roots in this case by transforming the objective function into a polynomial:

using PolynomialRoots
roots([-1.5e-6,-0.3,-15000,0,1])

produces

4-element Array{Complex{Float64},1}:
     122.47449713915809 - 0.0im
    -122.47447713915808 + 0.0im
 -1.0000000813048448e-5 + 0.0im
  -9.999999186951818e-6 + 0.0im

I would love a way to identify the pair of roots around the pole at x = -1e-5 without knowing the exact form of the objective function.

EDIT2: Trying out Roots.jl :

using Roots
f(x) = -15000 + x^4 / (x+1e-5)^2

find_zero(f,0.0) # finds +122... root
find_zero(f,(-1e-4,0.0)) # error, not a bracketing interval
find_zeros(f,-1e-4,0.0) # finds 0-element Array{Float64,1}
find_zeros(f,-1e-4,0.0,no_pts=6) # finds root slightly less than -1e-5
find_zeros(f,-1e-4,0.0,no_pts=10) # finds 0-element Array{Float64,1}, sensitive to value of no_pts

I can get find_zeros to work, but it's very sensitive to the no_pts argument and the exact values of the endpoints I pick. Doing a loop over no_pts and taking the first non-empty result might work, but something more deterministic to converge would be preferable.

EDIT3 : Here's applying the tanh transformation suggested by Bogumił

using NLsolve
function f_tanh!(F,x)
    x = x[1]
    x = -1e-4 * (tanh(x)+1) / 2
    F[1] = -15000 + x^4 / (x+1e-5)^2
end

nlsolve(f_tanh!,[100.0]) # doesn't converge
nlsolve(f_tanh!,[1e5]) # doesn't converge

using Roots
function f_tanh(x)
    x = -1e-4 * (tanh(x)+1) / 2
    return -15000 + x^4 / (x+1e-5)^2
end

find_zeros(f_tanh,-1e10,1e10) # 0-element Array
find_zeros(f_tanh,-1e3,1e3,no_pts=100) # 0-element Array
find_zero(f_tanh,0.0) # convergence failed
find_zero(f_tanh,0.0,max_evals=1_000_000,maxfnevals=1_000_000) # convergence failed

EDIT4 : This combination of techniques identifies at least one root somewhere around 95% of the time, which is good enough for me.

using Peaks
using Primes
using Roots

# randomize pole location
a = 1e-4*rand()
f(x) = -15000 + x^4 / (x+a)^2

# do an initial sample to find the pole location
l = 1000
minval = -1e-4
maxval = 0
m = []
sample_r = []
while l < 1e6
    sample_r = range(minval,maxval,length=l)
    rough_sample = f.(sample_r)
    m = maxima(rough_sample)
    if length(m) > 0
        break
    else
        l *= 10
    end
end
guess = sample_r[m[1]]

# functions to compress the range around the estimated pole
cube(x) = (x-guess)^3 + guess 
uncube(x) = cbrt(x-guess) + guess
f_cube(x) = f(cube(x))

shift = l ÷ 1000
low = sample_r[m[1]-shift]
high = sample_r[m[1]+shift]
# search only over prime no_pts, so no samplings divide into each other
# possibly not necessary?
for i in primes(500)
    z = find_zeros(f_cube,uncube(low),uncube(high),no_pts=i)
    if length(z)>0
        println(i)
        println(cube.(z))
        break
    end
end

Answer 1

More comment could be given if you provided more information on your problem.

However in general:

1. It seems that your problem is univariate, in which case you can use Roots.jl where find_zero and find_zeros give the interface you ask for (i.e. allowing to specify the search region)
2. If a problem is multivariate you have several options how to do it in the problem specification for nlsolve (as it by default does not allow to specify a bounding box AFAICT). The simplest is to use variable transformation. E.g. you can apply a ai * tanh(xi) + bi transformation selecting ai and bi for each variable so that it is bounded to the desired interval

The first problem you have in your definition is that the way you define f it never crosses 0 near the two roots you are looking for because Float64 does not have enough precision when you write 1e-5. You need to use greater precision of computations:

julia> using Roots

julia> f(x) = -15000 + x^4 / (x+1/big(10.0^5))^2
f (generic function with 1 method)

julia> find_zeros(f,big(-2*10^-5), big(-8*10^-6), no_pts=100)
2-element Array{BigFloat,1}:
 -1.000000081649671426108658262468117284940444265467160592853348997523986352593615e-05
 -9.999999183503552405580084054429938261707450678661727461293670518591720605751116e-06

and set no_pts to be sufficiently large to find intervals bracketing the roots.