有没有办法限制 NLsolve 在 Julia 中搜索的区域?
Is there any way to bound the region searched by NLsolve in Julia?
我试图找到一个非线性(大致为四次)方程的根。
该方程始终有四个根,其中一对接近于零,一个大的正根和一个大的负根。我想确定任何一个接近零的根,但是 nlsolve
,即使初始猜测非常接近这些根,似乎也总是收敛于大的正根或负根。
函数图基本上看起来像一个恒定的负值,在零附近有一个(非常窄的)偶数极点,并逐渐上升到在大的正根和负根处过零。
有什么方法可以限制 nlsolve
搜索的区域,或者做些什么让它对我的函数中这个极点的存在更加敏感?
编辑:
下面是重现问题的一些示例代码:
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])
作为输出,我得到:
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
在这种情况下,我们可以通过将 objective 函数转换为多项式来找到确切的根:
using PolynomialRoots
roots([-1.5e-6,-0.3,-15000,0,1])
产生
4-element Array{Complex{Float64},1}:
122.47449713915809 - 0.0im
-122.47447713915808 + 0.0im
-1.0000000813048448e-5 + 0.0im
-9.999999186951818e-6 + 0.0im
我希望有一种方法可以在不知道 objective 函数的确切形式的情况下识别 x = -1e-5 处极点周围的一对根。
编辑2:
尝试 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
我可以让 find_zeros
工作,但它对 no_pts
参数和我选择的端点的确切值非常敏感。在 no_pts
上进行循环并获取第一个非空结果可能会起作用,但更可取的是更确定性的收敛。
编辑3:
这里应用了 Bogumił
建议的 tanh 变换
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:这种技术组合在大约 95% 的时间内至少识别出一个根,这对我来说已经足够了。
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
如果您提供更多关于您的问题的信息,我们会给出更多评论。
但总的来说:
- 看来你的问题是单变量的,在这种情况下你可以使用 Roots.jl 其中
find_zero
和 find_zeros
给出你要求的界面(即允许指定搜索区域)
- 如果问题是多变量的,那么在
nlsolve
的问题规范中,您有多种选择(因为默认情况下不允许指定边界框 AFAICT)。最简单的是使用变量转换。例如。您可以应用 ai * tanh(xi) + bi
转换,为每个变量选择 ai
和 bi
,以便将其限制在所需的区间
您在定义中遇到的第一个问题是您定义 f
的方式永远不会越过 0
您正在寻找的两个根附近,因为 Float64
没有足够的写 1e-5
时的精度。您需要使用更高的计算精度:
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
并将 no_pts
设置得足够大以找到包含根的区间。
我试图找到一个非线性(大致为四次)方程的根。
该方程始终有四个根,其中一对接近于零,一个大的正根和一个大的负根。我想确定任何一个接近零的根,但是 nlsolve
,即使初始猜测非常接近这些根,似乎也总是收敛于大的正根或负根。
函数图基本上看起来像一个恒定的负值,在零附近有一个(非常窄的)偶数极点,并逐渐上升到在大的正根和负根处过零。
有什么方法可以限制 nlsolve
搜索的区域,或者做些什么让它对我的函数中这个极点的存在更加敏感?
编辑: 下面是重现问题的一些示例代码:
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])
作为输出,我得到:
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
在这种情况下,我们可以通过将 objective 函数转换为多项式来找到确切的根:
using PolynomialRoots
roots([-1.5e-6,-0.3,-15000,0,1])
产生
4-element Array{Complex{Float64},1}:
122.47449713915809 - 0.0im
-122.47447713915808 + 0.0im
-1.0000000813048448e-5 + 0.0im
-9.999999186951818e-6 + 0.0im
我希望有一种方法可以在不知道 objective 函数的确切形式的情况下识别 x = -1e-5 处极点周围的一对根。
编辑2:
尝试 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
我可以让 find_zeros
工作,但它对 no_pts
参数和我选择的端点的确切值非常敏感。在 no_pts
上进行循环并获取第一个非空结果可能会起作用,但更可取的是更确定性的收敛。
编辑3: 这里应用了 Bogumił
建议的 tanh 变换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:这种技术组合在大约 95% 的时间内至少识别出一个根,这对我来说已经足够了。
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
如果您提供更多关于您的问题的信息,我们会给出更多评论。
但总的来说:
- 看来你的问题是单变量的,在这种情况下你可以使用 Roots.jl 其中
find_zero
和find_zeros
给出你要求的界面(即允许指定搜索区域) - 如果问题是多变量的,那么在
nlsolve
的问题规范中,您有多种选择(因为默认情况下不允许指定边界框 AFAICT)。最简单的是使用变量转换。例如。您可以应用ai * tanh(xi) + bi
转换,为每个变量选择ai
和bi
,以便将其限制在所需的区间
您在定义中遇到的第一个问题是您定义 f
的方式永远不会越过 0
您正在寻找的两个根附近,因为 Float64
没有足够的写 1e-5
时的精度。您需要使用更高的计算精度:
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
并将 no_pts
设置得足够大以找到包含根的区间。