NLsolve not able to render solution

General Usage
,
1

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. :slight_smile:

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

2

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
3

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