Complex root finder for a general function f(z)?

IntervalRootFinding.jl should be able to do this but needs some work in this direction. Please do file issues on the repo.

1 Like

Thanks, but nlsolve can only find one root. The initial guess determines which root is found.

Thank you, I could convert it to a bivariate equation as Tamas suggested. But ForwardDiff (for autodiff) does not accept complex variables. Maybe this is the point that needs further work?

Good point. For now you would have to do Bisection first.

Yes, that’s how local rootfinders work. Rewritten, your problem is not different from a bivariate one, all caveats and trade-offs apply. In particular, if you need all roots, consider the interval methods suggested by @dpsanders, they are very powerful and well worth the setup cost.

Although you mention a “general complex equation in a certain complex region”, I would still be interested if you are motivated by a certain subclass with more structure to it.

For example, it often happens that such problems originate as nonlinear eigenvalue problems and it is desired to find all roots in a vertical strip, maybe for (in)stability considerations of some underlying dynamical system. In that case, NEP-PACK, also discussed here, may be interesting.

(I cannot help you with its use, but @jarl or his co-authors mentioned in that thread may be able to.)

3 Likes

Yes. NEP-PACK may be an efficient option if your nonlinear equation is actually det(M(s))=0 where M(s) is a square matrix. It may be a convenient option also for other nonlinear scalar analytic functions.

Scalar problem solving: If your problem can be expressed as a short sum of products of analytic functions (preferrably with matrix function implementations native in Julia) you can do something like this

julia> using NonlinearEigenproblems;
julia> f1= λ-> one(λ);
julia> f2=s->s^2;
julia> f3=s->sin(s^2);
julia> nep=SPMF_NEP([ones(1,1),ones(1,1),ones(1,1)],[f1,f2,f3]);
julia> (λ,v)=iar(nep,neigs=3,maxit=200);
julia> λ
3-element Array{Complex{Float64},1}:
 -7.751963555148087e-16 + 0.7148240548474641im
 -1.474965612325051e-15 - 0.714824054847463im 
    -0.4784100559939008 + 2.0871298203195203im
julia> 1+λ[1]^2+sin(λ[1]^2)
-2.3869795029440866e-15 - 2.0749569815028643e-15im

In the scalar case iar just reduces to dynamically computing the roots of the truncated Taylor expansion (for non-scalar it becomes a Kyrlov method). It’s based on computing derivatives around zero via the matrix function implementation of f1,f2 and f3. Note that sin(A) where Ais a matrix will give you the matrix function not the elementwise product (I think this is an awesome Julia design decision btw) which can be used to compute derivatives. So it’s basically a convenient way of computing roots of truncated Taylor series.

If it cannot be expressed using functions expressed with matrix function implementations, you may still be able to compute roots using some other NEP-PACK solvers, e.g. the contour integral methods contour_beyn and contour_block_SS, but then you need to tune some parameters. See list of NEP solvers: https://nep-pack.github.io/NonlinearEigenproblems.jl/methods/

6 Likes

Thanks for the suggestion. I am concerning about the solution of an equation such as the following:

(1+a*cosφ+b*cos2φ)^2 + (1+c*sinφ+d*sin2φ)^2 +1 = 0

where a,b,c,d are real constants, and I need the solution of φ in the region bound by the real part [- π, π] and imaginary part [0,50]. It seems to me that NEP-PACK is not supposed to solve such equations. Correct me if I am wrong. Thanks for the suggestion anyway!

Thank you very much for your detailed instructions. However, it seems that my problem does not seem to fulfill det(M(s))=0. I am concerning about the solution of equations such as that in the above post. Is NEP-PACK applicable in such cases? Thanks!

If you have a scalar equation of a single complex variable, you can treat it as the determinant of a 1x1 matrix M and use a nonlinear eigensolver routine.

2 Likes

Your problem seems not so difficult with this solution:

julia> using NonlinearEigenproblems
julia> a=3; b=4; c=-2; d=1;
julia> f1=s->one(s);
julia> f2=s->(one(s)+a*cos(s)+b*cos(2*s))^2
julia> f3=s->(one(s)+c*sin(s)+d*sin(2*s))^2
julia> nep=SPMF_NEP([ones(1,1),ones(1,1),ones(1,1)],[f1,f2,f3],check_consistency=false);
julia> (λ,v)=iar(nep,neigs=3,maxit=200);
3-element Array{Complex{Float64},1}:
 -1.0445788612507017 + 0.21155338513375282im
 -1.0445788612507032 - 0.21155338513375396im
   1.096857687197039 - 0.10586475855207848im
julia> φ=λ[1];  # Let's check solution 1 out of the three computed solutions
julia> (1+a*cos(φ)+b*cos(2φ))^2 + (1+c*sin(φ)+d*sin(2φ))^2 + 1
-2.682298827494378e-13 + 6.106226635438361e-14im

I would suggest a contour integral method (with ellipse) for this region, or provide a target for iar.

3 Likes

I see. Thanks!

Thank you. I appreciate your help. There are however two more questions:

  1. Why do you usually use three terms? If I just use 1 term, the calculation becomes much slower. Why? See the code below:
julia> using NonlinearEigenproblems
julia> a=3; b=4; c=-2; d=1;
julia> f(s) = one(s) + (one(s)+a*cos(s)+b*cos(2*s))^2 + (one(s)+c*sin(s)+d*sin(2*s))^2;
julia> nep=SPMF_NEP([ones(1,1)],[f],check_consistency=false);
julia> (λ,v)=iar(nep,neigs=3,maxit=200)
ERROR: No convergence: 'Number of iterations exceeded. maxit=200.Try to change the inner_solver_method for better performance.' eigenvalue approx:Complex{Float64}[-9.22885104278323 - 4.72099864915099im, -10.82326593488565 + 5.283169733939192im, -12.42691060068113 + 5.78799699741796im], errmeasure:[0.9999999999999998, 0.9999999999999998, 0.9999999999999998]
Stacktrace:
 [1] #iar#35(::Type{IterativeSolvers.DGKS}, ::Int64, ::FactorizeLinSolverCreator{Number,Any}, ::Float64, ::Int64, ::DefaultErrmeasure{StandardSPMFErrmeasure{Float64}}, ::Complex{Float64}, ::Complex{Float64}, ::Array{Float64,1}, ::Int64, ::Int64, ::Bool, ::DefaultInnerSolver, ::Int64, ::typeof(iar), ::Type{Complex{Float64}}, ::SPMF_NEP{Array{Float64,2},Complex{Float64}}) at /home/nih/.julia/packages/NonlinearEigenproblems/1hzkn/src/method_iar.jl:177
 [2] (::NonlinearEigenproblems.NEPSolver.var"#kw##iar")(::NamedTuple{(:neigs, :maxit),Tuple{Int64,Int64}}, ::typeof(iar), ::Type{Complex{Float64}}, ::SPMF_NEP{Array{Float64,2},Complex{Float64}}) at ./none:0
 [3] #iar#34 at /home/nih/.julia/packages/NonlinearEigenproblems/1hzkn/src/method_iar.jl:46 [inlined]
 [4] (::NonlinearEigenproblems.NEPSolver.var"#kw##iar")(::NamedTuple{(:neigs, :maxit),Tuple{Int64,Int64}}, ::typeof(iar), ::SPMF_NEP{Array{Float64,2},Complex{Float64}}) at ./none:0
 [5] top-level scope at REPL[67]:1

  1. I do not know how many roots there are a-priori for a general f(z). What should I do in that case?
    Thanks again!
  1. In this case the difference btw the iar application to my three terms NEP and your one term NEP is due to the default termination criteria. The two problems are equivalent, but their backward error are different, which is default error measure. Short fix (although I don’t see why you think this is a problem): if you add the kwarg errmeasure=ResidualErrmeasure(nep) to the iar call you should get more or less the same.

  2. If you set neigs=Inf and tol=Inf as kwarg in the iar-call, it will return whatever approximations are available after maxit iterations without throwing an error. You will however need to decrease maxit. (Since taking more Taylor terms does not always give better approximations due to round-off errors when computing roots of polynomials.) After that you will need to check which returned values are actually solutions:

julia> (λ,v)=iar(nep,neigs=Inf,tol=Inf,maxit=100);
julia> λ[abs.(f.(λ)).<1e-4]  # Candidate solutions 
8-element Array{Complex{Float64},1}:
 -1.0445788613522853 + 0.21155338512503322im
 -1.0445788613522855 - 0.2115533851250334im 
    1.09685768740707 + 0.10586475832144271im
    1.09685768740707 - 0.10586475832144267im
  2.6067446967142778 - 0.23816604552564838im
  -2.659023659901437 - 0.5456405751277179im 
   2.606744696714269 + 0.2381660455256412im 
 -2.6590236599014325 + 0.5456405751277291im 
1 Like

The only class of methods (I believe) which can guarantee (in the best case) that all the roots have been found are those based on interval arithmetic.

Using the method I outlined in https://github.com/JuliaDiffEq/ModelingToolkit.jl/issues/256, I found symbolic expressions for the real and complex parts of your function using ModelingToolkit.jl.

I used the same values of a etc. as in the posts above for comparison.
Then we just call the roots function from IntervalRootFinding.jl.
The option Krawczyk is one particular type of method for proving existence and uniqueness of roots using interval methods.

julia> using IntervalRootFinding, IntervalArithmetic, StaticArrays

julia> a=3; b=4; c=-2; d=1;

julia> f( (w, z) ) = SVector( (1 + (((1 + a * (cos(w) * cosh(z))) + b * (cos(2w) * cosh(2z))) * ((1 + a * (cos(w) * cosh(z))) + b * (cos(2w) * cosh(2z))) - (a * (-(sin(w)) * sinh(z)) + b * (-(sin(2w)) * sinh(2z))) * (a * (-(sin(w)) * sinh(z)) + b * (-(sin(2w)) * sinh(2z))))) + (((1 + c * (sin(w) * cosh(z))) + d * (sin(2w) * cosh(2z))) * ((1 + c * (sin(w) * cosh(z))) + d * (sin(2w) * cosh(2z))) - (c * (cos(w) * sinh(z)) + d * (cos(2w) * sinh(2z))) * (c * (cos(w) * sinh(z)) + d * (cos(2w) * sinh(2z)))), (((1 + a * (cos(w) * cosh(z))) + b * (cos(2w) * cosh(2z))) * (a * (-(sin(w)) * sinh(z)) + b * (-(sin(2w)) * sinh(2z))) + (a * (-(sin(w)) * sinh(z)) + b * (-(sin(2w)) * sinh(2z))) * ((1 + a * (cos(w) * cosh(z))) + b * (cos(2w) * cosh(2z)))) + (((1 + c * (sin(w) * cosh(z))) + d * (sin(2w) * cosh(2z))) * (c * (cos(w) * sinh(z)) + d * (cos(2w) * sinh(2z))) + (c * (cos(w) * sinh(z)) + d * (cos(2w) * sinh(2z))) * ((1 + c * (sin(w) * cosh(z))) + d * (sin(2w) * cosh(2z)))) )

julia> X = (-π..π) × (0..50)   # define the box as an IntervalBox
[-3.1416, 3.1416] × [0, 50]

julia> @time roots(f, X, Krawczyk, 1e-6)
  0.961258 seconds (4.58 M allocations: 209.121 MiB, 4.71% gc time)

4-element Array{Root{IntervalBox{2,Float64}},1}:
 Root([2.60674, 2.60675] × [0.238166, 0.238167], :unique)
 Root([-1.04458, -1.04457] × [0.211553, 0.211554], :unique)
 Root([1.09685, 1.09686] × [0.105864, 0.105865], :unique)
 Root([-2.65903, -2.65902] × [0.54564, 0.545642], :unique)

If all the code worked correctly, this should be a guarantee that there are precisely these 4 roots in the box! I’m not sure what happened to the other ones found in the post above. Presumably (hopefully) there are no actual roots there. EDIT: They have negative imaginary part and so are excluded from this calculation (as they were not requested in the OP). Doing the same search with imaginary part ranging in -50..50 successfully finds those extra 4 roots (which are complex conjugates of the ones found here).

Note that this is using the nonrecursive_powers branch of IntervalArithmetic.jl. If you use the released version you may (or may not) see a 10x or more slowdown.

4 Likes

Ah I see, the extra roots found in the previous post have negative imaginary part, which is excluded in my calculation.

By the way, you can get floating-point vectors instead of interval vectors as follows:

@time rts = roots(f, X, Krawczyk, 1e-6)

julia> mid.(interval.(rts))
4-element Array{SArray{Tuple{2},Float64,1,2},1}:
 [2.606744829958056, 0.23816603785182655]
 [-1.0445788627812245, 0.2115533847809708]
 [1.0968576871969855, 0.10586475855211779]
 [-2.6590236564036043, 0.5456406678966927]
1 Like

I generally use the well known Muller method, with deflation, and root polishing, for finding complex roots in Matlab. I have done a Google search, and can’t find the equivalent in Julia. The method is well described in “Numerical Recipes in C” (a Classic on numerical methods); it doesn’t involve a lot of coding; I wonder if someone could be enticed to put together a package. It really works well.

@kw_martin If you know the method well then maybe you would like to code up a version? I am definitely interested. How does deflation “know when to stop” for functions that are not polynomials?

1 Like

David, deflation is normally only used for polynomials, as far as
I know. I have personally only used Muller for finding polynomials
in filter design approximations where I knew how many roots to
look for. For this application, it worked well. For other
functions, if you have multiple roots, you probably need to know a
priori how many roots there are. For filter design approximation,
I am currently using Matlab where I already have a heavily
modified Muller routine; so it doesn’t make sense for me to code
it up in Julia (I’m guessing about a week of work - maybe more to
make it clean and a registered package which I have never done -
I’m new at Julia and struggling with a lot of issues like vectors,
real slow interpretation during debug, etc.). I see Julia as being
useful for realizing filters, but I don’t see it being used for
the design of filters.

1 Like