# 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[1] = kk[1] * x[1]
rate[2] = kk[2] * x[2]^2
rate[3] = kk[3] * x[3]

F[1] = 0.5*rate[3] - rate[1]
F[2] = rate[1] - rate[2]
F[3] = rate[1] + rate[2] - rate[3]
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[1] = kk[1] * x[1]
rate[2] = kk[2] * x[2]^2
rate[3] = kk[3] * x[3]

F[1] = 0.5*rate[3] - rate[1]
F[2] = rate[1] - rate[2]
F[3] = rate[1] + rate[2] - rate[3]
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[2] 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[1], x[2] and x[3] respectively, and there is another constraint that 2*x[1] + x[2] + x[3] = 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[3]*x[3]) - (kk[1]*x[1]) == 0
(kk[1]*x[1]) - (kk[2]*x[2]^2) == 0
(kk[1]*x[1]) + (kk[2]*x[2]^2) - (kk[3]*x[3]) == 0

which simplifies to:

2*kk[1]*x[1] - kk[3]*x[3] == 0
kk[1]*x[1] - kk[2]*x[2]^2 == 0
kk[1]*x[1] + kk[2]*x[2]^2 - kk[3]*x[3] == 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[3] = 2*x[1] + x[2] + x[3] - 1.0, I got the same solutions with what I got from DifferentialEquations.

2 Likes