non-linear solver always produces zero residual

Question

I am learning how to solve non-linear equations and Julia (1.3.1) in general and want to ask how we should use NLsolve.

As a first step, I try the following:

using NLsolve
uni = 1
z = 3
x = [3,2,4,5,6]
y = z .+ x
print(y)
function g!(F, x)
    F[1:uni+4] = y .- z - x
end
nlsolve(g!, [0.5,0.9,1,2,3])

And, I confirm it works as below:

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.5, 0.9, 1.0, 2.0, 3.0]
 * Zero: [2.999999999996542, 2.000000000003876, 4.000000000008193, 4.999999999990685, 5.999999999990221]
 * Inf-norm of residuals: 0.000000
 * Iterations: 2
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 3
 * Jacobian Calls (df/dx): 3

Then, I try a more complicated model as following

using SpecialFunctions, NLsolve, Random
Random.seed!(1234)
S = 2

# Setting parameters
ε = 1.3 
β = 0.4 
γ = gamma((ε-1)/ε) 
T = rand(S) 
E = rand(S)
B = rand(S)
w = rand(S)
Q = rand(S)
d = [1 2 ;2 1 ]

# Construct a model

rvector = T.*Q.^(1-β).*B.^ε
svector = E.* w.^ε


Φ_all = (sum(sum(rvector * svector' .* d )))
π = rvector * svector' .* d ./ Φ_all 

# These two are outcome the model 
πR = (sum(π,dims=1))' 
πM = sum(π,dims=2) 

# E is now set as unknown and we want to estimate it given the outcome of the model

function f!(Res, Unknown)


    rvector = T.*Q.^(1-β).*B.^ε
    svector = Unknown[1:S].* w.^ε 


    Φ_all = (sum(sum(rvector * svector' .* d )))
    π = rvector * svector' .* d ./ Φ_all
    Res = ones(S+S)
    Res[1:S] = πR - (sum(π,dims=1))' 
    Res[S+1:S+S] = πM - sum(π,dims=2) 

end


nlsolve(f!, [0.5,0.6])

and this code produces a weird result as follows.

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.5, 0.6]
 * Zero: [0.5, 0.6]
 * Inf-norm of residuals: 0.000000
 * Iterations: 0
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 1
 * Jacobian Calls (df/dx): 1

So, essentially, the function always yields 0 as the return, and thus the initial input always becomes the solution. I cannot understand why this does not work and also why this behaves differently from the first example. Could I have your suggestion to fix it?

Answer 1

You need to update Res in place. Now you do

Res = ones(S+S)

which shadows the input Res. Just update Res directly.