IntervalRootFinding.jl should be able to do this but needs some work in this direction. Please do file issues on the repo.
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 tradeoffs 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, NEPPACK, also discussed here, may be interesting.
(I cannot help you with its use, but @jarl or his coauthors mentioned in that thread may be able to.)
Yes. NEPPACK 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> λ
3element Array{Complex{Float64},1}:
7.751963555148087e16 + 0.7148240548474641im
1.474965612325051e15  0.714824054847463im
0.4784100559939008 + 2.0871298203195203im
julia> 1+λ[1]^2+sin(λ[1]^2)
2.3869795029440866e15  2.0749569815028643e15im
In the scalar case iar
just reduces to dynamically computing the roots of the truncated Taylor expansion (for nonscalar 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 A
is 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 NEPPACK 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://neppack.github.io/NonlinearEigenproblems.jl/methods/
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 NEPPACK 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 NEPPACK 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.
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);
3element 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.682298827494378e13 + 6.106226635438361e14im
I would suggest a contour integral method (with ellipse) for this region, or provide a target for iar
.
I see. Thanks!
Thank you. I appreciate your help. There are however two more questions:
 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] toplevel scope at REPL[67]:1
 I do not know how many roots there are apriori for a general f(z). What should I do in that case?
Thanks again!

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. 
If you set
neigs=Inf
andtol=Inf
as kwarg in the iarcall, it will return whatever approximations are available aftermaxit
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 roundoff 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.(λ)).<1e4] # Candidate solutions
8element 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
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, 1e6)
0.961258 seconds (4.58 M allocations: 209.121 MiB, 4.71% gc time)
4element 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.
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 floatingpoint vectors instead of interval vectors as follows:
@time rts = roots(f, X, Krawczyk, 1e6)
julia> mid.(interval.(rts))
4element Array{SArray{Tuple{2},Float64,1,2},1}:
[2.606744829958056, 0.23816603785182655]
[1.0445788627812245, 0.2115533847809708]
[1.0968576871969855, 0.10586475855211779]
[2.6590236564036043, 0.5456406678966927]
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?
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.