Root finding for (possibly) complex valued functions with Nonlinearsolve.jl

Hello.

I need to translate some matlab code into Julia, where i need to find the roots in X of the function:

function problem(X :: AbstractVector{T}, A, β, Ω, α, ε, θ, σ,L) where T

    N = length(α)
    p = @view X[1:N]
    y = @view X[N+1:end]
    
    Out :: Vector{T} = zeros(eltype(X),2*N)
    
    q = (Ω * p .^ (1-θ)) .^ (1 / (1 - θ))
    w = p .* (A .^ ((ε - 1)/ε)) .* (α .^ (1 / ε)) .* (y .^ (1/ε)) .* L .^ (-1/ε)
    C = w' * L
  
    Out[1:N] = p - (A .^ (ε - 1) .* (α .* w .^ (1- ε) + (1 .- α) .* q .^ (1 - ε))) .^ (1/(1-ε))
    Out[N+1:end] = y' - y' * diagm(p)^ε * diagm(A)^(ε-1) * diagm(q)^(θ-ε) * diagm(1 .- α) * Ω * diagm(p)^(-θ) - β'*diagm(p)^(-σ)*C
    
    return Out
end

β, Ω, α, ε, θ, σ and L are fixed parameters. A is a random variable.
Nonlinearsolve.jl an NLsolve.jl both are able to give exact solutions for norm A close enought to one, however for smaller A the problem function returns complex values, and here I run into problems.
If I allow complex values either ForwardDiff stops working, or I get a Innexact Error by the solver. I tried simply squaring Out, however I get very unexact results then. Fsolve works fine it matlab for the same problem.

Thank you in advance for your answers.

Have you tried NewtonRhapson(autodiff=false)?

2 Likes

Thank you, that did do the trick.