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

Specific Domains Optimization (Mathematical)
, ,
1

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.

2

Have you tried NewtonRhapson(autodiff=false)?

2 Likes
3

Thank you, that did do the trick.

4

When I try this I get an error. The following code

params = (; a, b, m, n, k, ηwn)
Δs0 = @SVector[0.0, 0.0, 0.0, 0.0]
prob = Snls.NonlinearProblem(WaveguideModes.nlsolvedkxadkyb, Δs0, params)
sol = Snls.solve(prob, Snls.NewtonRaphson(autodiff=false))

results in

ERROR: MethodError: no method matching select_jacobian_autodiff(::SciMLBase.NonlinearProblem{…}, ::Bool)
The function `select_jacobian_autodiff` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  select_jacobian_autodiff(::SciMLBase.AbstractNonlinearProblem, ::Nothing)
   @ NonlinearSolveBase C:\Users\peter\.julia\packages\NonlinearSolveBase\lDNHI\src\autodiff.jl:87
  select_jacobian_autodiff(::SciMLBase.AbstractNonlinearProblem, ::ADTypes.AbstractADType)
   @ NonlinearSolveBase C:\Users\peter\.julia\packages\NonlinearSolveBase\lDNHI\src\autodiff.jl:75

Stacktrace:
 [1] macro expansion
   @ C:\Users\peter\.julia\packages\Setfield\ZezIj\src\sugar.jl:198 [inlined]
 [2] __init(::SciMLBase.NonlinearProblem{…}, ::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, maxiters::Int64, abstol::Nothing, reltol::Nothing, maxtime::Nothing, termination_condition::Nothing, internalnorm::Function, linsolve_kwargs::@NamedTuple{}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{})
   @ NonlinearSolveFirstOrder C:\Users\peter\.julia\packages\NonlinearSolveFirstOrder\5esMl\src\solve.jl:130
 [3] __solve(::SciMLBase.NonlinearProblem{…}, ::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…}; kwargs::@Kwargs{…})
   @ NonlinearSolveBase C:\Users\peter\.julia\packages\NonlinearSolveBase\lDNHI\src\solve.jl:5
 [4] solve_call(_prob::SciMLBase.NonlinearProblem{…}, args::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
   @ DiffEqBase C:\Users\peter\.julia\packages\DiffEqBase\4uSBa\src\solve.jl:127
 [5] solve_up(prob::SciMLBase.NonlinearProblem{…}, sensealg::Nothing, u0::StaticArraysCore.SVector{…}, p::@NamedTuple{…}, args::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…}; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{…})
   @ DiffEqBase C:\Users\peter\.julia\packages\DiffEqBase\4uSBa\src\solve.jl:665
 [6] solve(prob::SciMLBase.NonlinearProblem{…}, args::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
   @ DiffEqBase C:\Users\peter\.julia\packages\DiffEqBase\4uSBa\src\solve.jl:638
 [7] solve(prob::SciMLBase.NonlinearProblem{…}, args::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…})
   @ DiffEqBase C:\Users\peter\.julia\packages\DiffEqBase\4uSBa\src\solve.jl:604
 [8] top-level scope
   @ c:\Users\peter\Documents\julia\packages\WaveguideModes\test\junk.jl:17
Some type information was truncated. Use `show(err)` to see complete types.

Julia version and runtime environment:

julia> versioninfo()
Julia Version 1.11.6
Commit 9615af0f26 (2025-07-09 12:58 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 28 × Intel(R) Core(TM) i7-14700
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 28 default, 0 interactive, 14 GC (on 28 virtual cores)

Package environment:

(WaveguideModes) pkg> st
Project WaveguideModes v0.1.0
Status `C:\Users\peter\Documents\julia\packages\WaveguideModes\Project.toml`
  [6e696c72] AbstractPlutoDingetjes v1.3.2
  [7d9f7c33] Accessors v0.1.42
  [b21f74c0] FunctionZeros v1.0.0 `https://github.com/simonp0420/FunctionZeros.jl#`
  [d1286fa9] MetalSurfaceImpedance v1.0.1
  [8913a72c] NonlinearSolve v4.10.0
  [5ad8b20f] PhysicalConstants v0.2.4
⌅ [08abe8d2] PrettyTables v2.4.0
  [90137ffa] StaticArrays v1.9.15
  [1986cc42] Unitful v1.24.0
  [de0858da] Printf v1.11.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`
5

Sorry that’s the old syntax. autodiff=AutoFiniteDiff() from ADTypes.jl

1 Like
6

:+1: Thanks for the rapid reply!