NLsolve always gives the initial guess as the final answer

Specific Domains Modelling & Simulations
1

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 ?

2

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

3

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)?

4

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
6

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