# NLsolve not able to render solution

Hi everyone,

I am trying to use `NLsolve` to solve a nonlinear equation and trying to compare it with `scipy.optimize.root`. However, as usual, the solver isn’t rendering the solution and my guess is that I am not using it properly. If you think that a similar question has already been answered, then feel free to point me to it.

I am passing the same initial guess in both `scipy.optimize.root` and `nlsolve`, so I do expect that the latter should converge.

### Julia code

``````using Optim, Plots, LaTeXStrings, NLsolve
gr()

I1 = (λ1, λ2) -> @. λ1^2 + λ2^2 + 1/(λ2^2 * λ1^2 )
I2 = (λ1, λ2) -> @. 0.5 * (I1(λ1, λ2)^2 - (λ1^4 + λ2^4 + 1/(λ2^4 * λ1^4 )))

cond1 = (λ1, λ2) -> @. I2(λ1, λ2)/I1(λ1, λ2) - 0.762 * I1(λ1, λ2) + 9.22

# solving condition 1
l2vals = Array(3.5^(-0.5):0.1:2)
function func!(F, l1)
F = cond1(l1, l2vals)
end

sol = nlsolve(func!, 3.5*ones(size(l2vals)), autodiff = :forward, method=:newton)
``````

which gives the following output

``````Results of Nonlinear Solver Algorithm
* Algorithm: Newton with line-search
* Starting Point: [3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5]
* Zero: [-Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf]
* Inf-norm of residuals: 3.500000
* Iterations: 1000
* Convergence: false
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: false
* Function Calls (f): 2
* Jacobian Calls (df/dx): 2
``````

### python

``````import numpy as np
from scipy.optimize import root

def I1(l1, l2):
return l1**2. + l2**2. + 1./(l1**2. * l2**2.)

def I2(l1, l2):
return 1./2 * (I1(l1, l2)**2. - l1**4. - l2 ** 4. - 1./(l1**4. * l2**4.))

def cond1(l1, l2):
I1v = I1(l1, l2)
I2v = I1(l1, l2)
res = I2v/I1v - 0.762*I1v + 9.22
return res

l2vals = np.linspace(3.5**(-0.5), 2., 100)
func = lambda x: cond1(x, l2vals)
sol = root(func, 3.5*np.ones_like(l2vals))
print(sol.x, "\n", sol.message, sol.fun)
``````

which works fine…

``````[3.58525758 3.58503458 3.58458982 3.58393993 3.5830994  3.58208088
3.58089545 3.57955286 3.5780617  3.57642958 3.57466324 3.57276867
3.5707512  3.56861559 3.56636612 3.56400659 3.56154043 3.55897072....
The solution converged.
[-2.14010143e-10 -2.06224371e-10 -1.95967687e-10 -1.83785431e-10
-1.70196301e-10 -1.55690572e-10 -1.40705225e-10 -1.25645272e-10
-1.10837561e-10 -9.65680869e-11 -8.30642222e-11 -7.04982739e-11
-5.89910343e-11 -4.86188867e-11 -3.94191346e-11 -3.13935544e-11 ....
``````

Could anyone point me to a solution?

Cheers

Ah, it was function.

``````function func!(F, l1)
@. F = cond1(l1, l2vals)
end
``````

or more precisely

``````function func!(F, l1)
F .= cond1(l1, l2vals)
end
``````

Since `I1, I2 and cond1` is all element-wise just make them scalar. At the end, broadcasted `setindex!` can be non-allocating. Like

``````I1 = (λ1, λ2) -> λ1^2 + λ2^2 + 1/(λ2^2 * λ1^2 )
I2 = (λ1, λ2) -> 0.5 * (I1(λ1, λ2)^2 - (λ1^4 + λ2^4 + 1/(λ2^4 * λ1^4 )))

cond1 = (λ1, λ2) -> I2(λ1, λ2)/I1(λ1, λ2) - 0.762 * I1(λ1, λ2) + 9.22

# solving condition 1
l2vals = Array(3.5^(-0.5):0.1:2)
function func!(F, l1)
F .= cond1.(l1, l2vals)  # LHS does not allocate
end
``````

But there are other allocation points like `I1, I2, cond1 and l2vals` all allocate as they are accessed because they are nonconstant global variables.

1 Like