 # NLsolve always gives the initial guess as the final answer

Hi there! I’m trying to solve a nonlinear system of equations using the `NLsolve` package. However, I kept getting the same answer as my initial guess as if the solver did not do any iteration at all. I did not use Jacobian as in my real applications an analytic Jobabian is not available. An example code:

``````kk = [2e-17, 8e-12, 3e-17]
initial_x = [0.1, 0.0, 0.8]

using NLsolve

function f!(F, x)
rate = zeros(3)
rate = kk * x
rate = kk * x^2
rate = kk * x

F = 0.5*rate - rate
F = rate - rate
F = rate + rate - rate
end

sol = nlsolve(f!, initial_x, ftol=1e-20)
``````

with the following outputs:

Results of Nonlinear Solver Algorithm

• Algorithm: Trust-region with dogleg and autoscaling
• Starting Point: [0.5, 0.0, 0.0]
• Zero: [0.5, 0.0, 0.0]
• Inf-norm of residuals: 0.000000
• Iterations: 1000
• Convergence: false
• |x - x’| < 0.0e+00: false
• |f(x)| < 1.0e-20: false
• Function Calls (f): 1
• Jacobian Calls (df/dx): 1

It seems that there has been 1000 iterations but somehow the solution (`sol.zero`) is still the same as the initial guess. I played around with the `ftol` parameter but that did not help.

On the other hand, if I used `DifferentialEquations` and integrated towards a steady-state solution, I got the solution I expected analytically:

``````using DifferentialEquations

function fode!(F, x, p, t)
rate = zeros(3)
rate = kk * x
rate = kk * x^2
rate = kk * x

F = 0.5*rate - rate
F = rate - rate
F = rate + rate - rate
end

tend = 1e17
tspan = (0.0, tend)
prob = ODEProblem(fode!, initial_x, tspan)
sol = solve(prob);
#sol = solve(prob, reltol=1e-7, abstol=1e-7);
xf = sol.u[end]
``````

What did I do wrong with `NLsolve `?

whats the correct answer?, the inicial vector `[0.1, 0.0, 0.8]` has any restrictions?

The correct answer is roughly xf ~ [0.3, 0, 0.4] where xf is a small number. With `DifferentialEquations`, I got [0.2997403048544782, 0.0008656504849742739, 0.3996537398060098] which is pretty close. The element of the vector must be in the range of [0, 0.5], [0, 1] and [0, 1] for x, x and x respectively, and there is another constraint that 2*x + x + x = 1.0, which should remain true if it is true initially. Is there any way to add constraints to the solver (e.g. specifying the allowed range of the solution)?

Are these equations independent?

You are trying to solve the system:

``````0.5*(kk*x) - (kk*x) == 0
(kk*x) - (kk*x^2) == 0
(kk*x) + (kk*x^2) - (kk*x) == 0
``````

which simplifies to:

``````2*kk*x - kk*x == 0
kk*x - kk*x^2 == 0
kk*x + kk*x^2 - kk*x == 0
``````

equation 1 minus equation 2 yields equation 3.

2 Likes

Thanks for the remark! As you pointed out, the three equations are not independent and so the Jacobian is singular, which is probably why `NLsolve` didn’t work. Once I replaced the third equation with `F = 2*x + x + x - 1.0`, I got the same solutions with what I got from `DifferentialEquations`.

2 Likes