Basic usage of NLsolve for scalar problems

Specific Domains Optimization (Mathematical)
1

I have troubles to use NLsolve package for solving the following simple scalar example:

0 = cos(x)

for x initialized to, say, 1.5 (the anticipated solution is thus pi/2, that is, something like 1.57). The code is

julia> using NLsolve

julia> function f!(r,x)
             r = cos.(x) # x (and r) are expected to be vectors, hence the broadcasting
       end
f! (generic function with 1 method)

julia> function j!(J,x)
             J = -sin.(x)
       end
j! (generic function with 1 method)

julia> sol = nlsolve(f!,j!,[1.5])
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [1.5]
 * Zero: [1.5]
 * 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

julia> sol.zero
1-element Array{Float64,1}:
 1.5

Not even a single iteration is run. What am I doing wrong? Thanks in advance.

2

maybe this?

function f!(r,x)
    r .= cos.(x) # x (and r) are expected to be vectors, hence the broadcasting
end

function j!(J,x)
    (s1,s2) = size(J) #the jacobian is a square matrix, not a vector. the analytical derivative is Diagonal(-sin(x))
    J .= zeros(s1,s1)
    for i in 1:s1
        J[i,i] = -sin(x[i])
    end
end
1 Like
3

Yup. The r = cos.(x) doesn’t actually change the value of r, it just reassigns it. For scalar problems, use Roots, though.

1 Like
4

Thanks to both of you. Indeed, the problem was in the difference between = and .=. The modified code

using NLsolve

function f!(r,x)
      r .= cos.(x) # x (and r) are expected to be vectors, hence the broadcasting
end

function j!(J,x)
      J .= -sin.(x)
end

sol = nlsolve(f!,j!,[1.5])

sol.zero

is functional.

I think I need to read a bit more about the .= stuff. I guess that this paragraph

Finally, the maximum efficiency is typically achieved when the output array of a vectorized operation is pre-allocated , so that repeated calls do not allocate new arrays over and over again for the results (see Pre-allocating outputs). A convenient syntax for this is X .= ... , which is equivalent to broadcast!(identity, X, ...) except that, as above, the broadcast! loop is fused with any nested “dot” calls. For example, X .= sin.(Y) is equivalent to broadcast!(sin, X, Y) , overwriting X with sin.(Y) in-place. If the left-hand side is an array-indexing expression, e.g. X[2:end] .= sin.(Y) , then it translates to broadcast! on a view , e.g. broadcast!(sin, view(X, 2:lastindex(X)), Y) , so that the left-hand side is updated in-place.

from the section " Dot Syntax for Vectorizing Functions" in the manual contains the explanation. I might need some quiet time to digest it.

And, of course, in this particular case, using Roots.jl is certainly the way to solve the problem but I wanted to understand this issue to learn Julia itself better. Thanks.

1 Like