I think I found a bug in the NonlinearSolve.jl package. Consider the problem
using NonlinearSolve, GLMakie
f(x, p) = tanh(x / 1e-3) + x - 0.5
x0 = 1.0
prob = NonlinearProblem(f, x0)
sol = solve(prob, NewtonRaphson(), abstol=1e-3)
which returns
julia> sol.u
1.5
julia> f(sol.u, nothing)
2.0
This is clearly not a converged solution, since it is far from a zero of the function. However, no error is thrown by the solver even though the Newton iterations are not converging. This non-convergence is due to the shape of this nonlinear function
in the case where the zero is in the steep part of the graph. For an initial condition on the regions of the function dominated by x, you just end up iterating between two different values and do not achieve convergence.
I understand that this is an expected behavior of the Newton Raphson method for this function, but I think NonlinearSolve.jl should throw and error instead of letting you think that convergence was achieved.
and it seems that the maximum of 1000 iterations has been reached (I assumed this number is the one by default on the package.
I was expecting an error similar to how it is handled when solving and ODE with OrdinaryDiffEq, where it is explicitly stated when the integration is not converging within tolerances for adaptive timesteps.
It stopped at max iterations and exited early as a failure. That’s what you’re seeing?
OrdinaryDiffEq also returns a retcode that describes the solution, the same retcode actually. Indeed we should make NonlinearSolve.jl throw a warning for this case, but that’s a feature request for the repo (or a PR implementation) and not really a Discourse discussion.
That would be nice, if I have time next week I may try to implement it myself. The problem is that even though the no convergence information was in the retcode, it did not seem obvious to me and I assumed that the nonlinear solver succeeded. It was not until I plotted the results that I noticed that truly there was no convergence, which I think it could be misleading when using the package.
Even better than a warning, I think that if the full description of the solution object was printed in the REPL would already be okay to see at a glance if convergence happened. Now I was only getting sol.u printed in the REPL.
