How to deal with vectors in NLsolve.jl

[khaledharizb]

Hi guys!

I didn’t find a way to deal with vector using pkg NLSolve in the following problem.

Let q=[x,y] and p=[u,v] and defining the functions f and g as f(q) = q / (√(q'*q))^3,\ g(p) = p. Now, giving q_0=[0.4,0] and p_0=[0,0.5]. I need to solve by Newton-Raphson method the following

F(q_0,p_0,q,p) = [q - p_0 + 0.5 \cdot f((q_0+q)/2) , p - q_0 - 0.5 \cdot g((p_0+p)/2)]

using ForwardDiff

f(q) = q / (√(q'*q))^3;
g(p) = p;   
    
F(q₀,p₀,q,p) = [q - p₀ + 0.5 * f((q₀+q)/2) , p - q₀ - 0.5 * g((p₀+p)/2)]
ini_gauss = [q₀,p₀]
nlsolve((q,p)-> F(q₀,p₀,q,p), ini_gauss, autodiff = :forward, method=:newton, ftol=1e-14) 

[Mason]

Your parameters can’t be a vector of vectors, you need to do something like

ini_gauss = [q₀..., p₀...]

i.e. you need to flatten this structure. You’ll also need to adjust that closure, i.e.

(q,p)-> F(q₀,p₀,q,p)

should be

qp -> F(q₀,p₀, view(qp, 1:2), view(qp, 3:4))

and similarly, F should return a flat vector instead of a vector of vectors:

function F(q₀,p₀,q,p) 
    p′ = q - p₀ + 0.5 * f((q₀+q)/2) 
    q′ = p - q₀ - 0.5 * g((p₀+p)/2)
    [p′..., q′...]
end

Note, the above code is not very optimized and there’s a lot that could be done to make it faster. However, at least with the initial conditions you gave, it seems this system of equations does not converge to a solution with any of the solvers I tried, so I wonder if there’s a problem with your initial conditions or your equations.

[khaledharizb]

It works, great! Many thanks @Mason