Not finding a zero for a fairly simple problem

[elbert5770]

This doesn't find the zero at u = 0.070720181679945, using either a bounded or unbounded method. Is this a problem with the setup or is it more fundamental?

using NonlinearSolve
using Plots

f(u) = 0.5/1.5*log.(u./(1.0.-u)) .- 2.0*u .+1.0

uspan = (0.02, 0.1) # brackets
probB = IntervalNonlinearProblem(f, uspan)
sol = solve(probB, Falsi())
@show sol.u

u0 = 0.06
p = 2.0
probN = NonlinearProblem(f, u0, p)
solver = solve(probN, NewtonRaphson(), reltol = 1e-9)
@show solver.u

@show f(0.070720181679945)

c = collect(0.01:0.001:0.99)
plot(c,f(c))

And the output is:

sol.u = 0.10000000000000002
solver.u = 0.5
f(0.070720181679945) = 6.661338147750939e-16

Here's MATLAB code that gives the correct answer using both a bounded and unbounded method:

function test
u = linspace(0.01,0.99);
for i = 1:length(u)
   y(i) = f(u(i));
end
%plot(u,y)
x0 = [0.02 0.1]
z1 = fzero(@f,x0)
z2 = fzero(@f,0.06)

function y = f(u)
y = 0.5/1.5*log(u./(1.0 - u)) - 2.0 .*u + 1.0;

[sjdaines]

I can't reproduce this (with NonlinearSolve v1.8.0) and I think there are two problems here + a style point with use of broadcasting.

1. A problem with the setup. f needs to be f(u, p), where the parameters p are unused in this case.

From a fresh REPL, you should get an error explaining this, which starts:

julia> probB = IntervalNonlinearProblem(f, uspan)
ERROR: All methods for the model function `f` had too few arguments. For example,
an ODEProblem `f` must define either `f(u,p,t)` or `f(du,u,p,t)`. This error
can be thrown if you define an ODE model for example as `f(u,t)`. The parameters
`p` are not optional in the definition of `f`! For more information on the required
number of arguments for the function you were defining, consult the documentation
for the `SciMLProblem` or `SciMLFunction` type that was being constructed.

...

Plausibly what has happened here is that you have previously defined a method for f earlier in the same REPL session which did have the correct signature (with two arguments), and it is this that is being called instead of the single-argument method.

2. The function has multiple roots, and the Newton solution in your example is finding a different root at u = 0.5 (although still not calling the correct f, so this is really a coincidence) which is adding to the confusion.

3. In Julia it is clearer to define f(u, p) as a scalar function of u, and then vectorize where needed (this doesn't look like it affects the solution here)

With a modified example, I get:

using NonlinearSolve
using Plots

f(u, p) = 0.5/1.5*log(u/(1.0-u)) - 2.0*u +1.0 

uspan = (0.02, 0.1) # brackets
probB = IntervalNonlinearProblem(f, uspan)
sol = solve(probB, Falsi())

u0 = 0.06
p = 2.0
probN = NonlinearProblem(f, u0, p)
solver = solve(probN, NewtonRaphson(), reltol = 1e-9)

f(0.070720181679945)

c = collect(0.01:0.001:0.99)
plot(c, f.(c, nothing))   # broadcast scalar f

where the relevant parts of the output are:

julia> sol = solve(probB, Falsi())
u: 0.07072018167994475

julia> solver = solve(probN, NewtonRaphson(), reltol = 1e-9)
u: 0.07072017898982122

[ChrisRackauckas]

Indeed, NonlinearSolve.jl never supported f(u) as the syntax, so it's unclear what the OP even ran but it's definitely not the function that's posted in the issue. In any case, the current default solvers solve this problem just fine:

using NonlinearSolve

f(u, p) = 0.5/1.5*log.(u./(1.0.-u)) .- 2.0*u .+1.0

uspan = (0.02, 0.1) # brackets
prob = IntervalNonlinearProblem(f, uspan)
sol = solve(prob) # 0.07072018167994477
@show sol.resid # 0.0

u0 = 0.06
p = 2.0
prob = NonlinearProblem(f, u0, p)
solver = solve(prob)
@show solver.u # 0.07072018167994477
@show sol.resid # 0.0

so closing.