Hello, I am new to Julia. I wonder if there is a Julia package that can obtain the complex root for a general complex equation f(z)=0? I am aware there are packages such as Polynomials and PolynomialRoots, but they can only find roots for polynomial functions. I am looking for something that can obtain all complex roots to a general complex equation f(z)=0 in a certain complex region. Any suggestions? Thanks a lot!
2 Likes
Solving a bivariate
x -> (y = f(Complex(x[1], x[2])); [real(y), imag(y)])
could work, eg with
2 Likes
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.
1 Like
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: NEP-Solvers · NEP-PACK
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:
- 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
- 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!
-
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=Infandtol=Infas kwarg in the iar-call, it will return whatever approximations are available aftermaxititerations 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 Complex functions · Issue #34 · JuliaSymbolics/Symbolics.jl · GitHub, 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.