# Complex root finder for a general function f(z)?
**URL:** https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989
**Category:** New to Julia
**Tags:** question, package, roots
**Created:** [March 14, 2020, 9:12pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989 "2020-03-14T21:12:26Z")
**Posts on this page:** 20
**Page:** 2
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 23, 2020, 4:38am UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/21 "2020-03-23T04:38:11Z")
@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?
---
### Author:  [@kw\_martin](https://discourse.julialang.org/u/kw_martin)
#### Post date: [March 23, 2020, 2:38pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/22 "2020-03-23T14:38:56Z")
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.
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 23, 2020, 6:33pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/23 "2020-03-23T18:33:17Z")
Well I’m mainly suggesting that you start by coding up the main algorithm and posting it here. That’s the quickest way to get useful feedback and improve iteratively, and that’s the fastest way to get a working implementation of the algorithm in Julia! It will also be a good learning experience.
It doesn’t need to be performant to start with – that will come later. And since it’s Julia, it will be easy to make it performant once the first implementation is available.
---
### Author:  [@hongchengni](https://discourse.julialang.org/u/hongchengni)
#### Post date: [March 25, 2020, 11:08am UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/24 "2020-03-25T11:08:58Z")
Thanks for your help!
---
### Author:  [@hongchengni](https://discourse.julialang.org/u/hongchengni)
#### Post date: [March 25, 2020, 11:19am UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/25 "2020-03-25T11:19:47Z")
Thanks for your help! That’s very useful. The only problem is that I need to manually get the analytical real and imaginary part beforehand. I have also used SymPy to get the real and imaginary parts:
```julia
using SymPy
a = symbols("a",real=true)
b = symbols("b",real=true)
c = symbols("c",real=true)
d = symbols("d",real=true)
φr = symbols("φr",real=true)
φi = symbols("φi",real=true)
expr = (1+a*cos(φr+1im*φi)+b*cos(2φr+2im*φi))^2 + (1+c*sin(φr+1im*φi)+d*sin(2φr+2im*φi))^2 + 1
exprexpand = expr.expand(complex=true)
println(real(exprexpand))
println(imag(exprexpand))
```
Unfortunately, the actual problem I’m trying to solve now is extremely complicated, and the complex expansion was running like for ever. Even if it gets out finally, the analytical expression will be extremely long, and its evaluation will take a long time. It would still be great if IntervalRootFinding.jl natively supports complex root finding.
---
### Author:  [@abulak](https://discourse.julialang.org/u/abulak)
#### Post date: [March 25, 2020, 12:30pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/26 "2020-03-25T12:30:06Z")
looking at [http://www.grad.hr/nastava/gs/prg/NumericalRecipesinC.pdf](http://www.grad.hr/nastava/gs/prg/NumericalRecipesinC.pdf) (p. 371) the Muller method seems to be
```julia
qq(xᵢ₋₂, xᵢ₋₁, xᵢ) = (xᵢ - xᵢ₋₁)/(xᵢ₋₁ - xᵢ₋₂)
function muller(P, lo, up, atol=eps(one(lo)))
xᵢ₋₂, xᵢ₋₁, xᵢ = lo, (lo+up)/2, up
# @show abs(xᵢ₋₂ - xᵢ)
while abs(xᵢ₋₂ - xᵢ) > atol
q = qq(xᵢ₋₂, xᵢ₋₁, xᵢ)
q² = q^2
q1 = q+one(q)
pᵢ = P(xᵢ)
pᵢ₋₁= P(xᵢ₋₁)
pᵢ₋₂= P(xᵢ₋₂)
A = q*pᵢ - q*q1*pᵢ₋₁ + q²*pᵢ₋₂
B = (q1+q)*pᵢ - q1^2*pᵢ₋₁ + q²*pᵢ₋₂
C = q1*pᵢ
den = let
Δ = sqrt(B^2 - 4A*C)
d⁺ = B + Δ
d⁻ = B - Δ
abs(d⁺) > abs(d⁻) ? d⁺ : d⁻
end
x = xᵢ - (xᵢ - xᵢ₋₁)*(2C/den)
xᵢ₋₂, xᵢ₋₁, xᵢ = xᵢ₋₁, xᵢ, x
# @show xᵢ₋₂, xᵢ₋₁, xᵢ
end
return xᵢ₋₁
end
f(x) = x^2-2
muller(f, 1.0, 2.0) # converges in 3 interations
@btime muller($f, 1.0, 2.0) # 71.284 ns (0 allocations: 0 bytes)
```
In which direction would you like to take it?
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 25, 2020, 2:31pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/27 "2020-03-25T14:31:38Z")
@abulak Did you translate their code? Numerical recipes has a terribly restrictive software license and it’s basically not allowed to use their code.
---
### Author:  [@abulak](https://discourse.julialang.org/u/abulak)
#### Post date: [March 25, 2020, 2:33pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/28 "2020-03-25T14:33:15Z")
@dpsanders there is no code for this method in the book quoted; this is literal translation of formulas
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 25, 2020, 3:30pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/29 "2020-03-25T15:30:02Z")
OK great, thanks!
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 25, 2020, 3:32pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/30 "2020-03-25T15:32:01Z")
Very nice code, by the way 😉
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 25, 2020, 3:33pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/31 "2020-03-25T15:33:57Z")
@abulak I’m getting `qq` not defined.
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 25, 2020, 5:32pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/32 "2020-03-25T17:32:39Z")
@hongchengni I have just realised that it should be possible to directly use the complex function as long as you provide the complex derivative, since from that we can extract the Jacobian of the real functions (real and imaginary parts). This may require some rewriting / refactoring of the package:
[https://github.com/JuliaIntervals/IntervalRootFinding.jl/issues/147](https://github.com/JuliaIntervals/IntervalRootFinding.jl/issues/147)
Unfortunately I believe that automatic differentiation with complex functions is still a grey area?
I believe that there is also an equivalent of the interval Newton method for working with complex intervals, but I am having trouble finding this in the literature. (I tried the naive thing with polynomials and it seemed to work.)
---
### Author:  [@abulak](https://discourse.julialang.org/u/abulak)
#### Post date: [March 25, 2020, 5:59pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/33 "2020-03-25T17:59:08Z")
sorry, qq got pasted on the 0-th line; i corrected now
---
### Author:  [@abulak](https://discourse.julialang.org/u/abulak)
#### Post date: [March 25, 2020, 9:05pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/34 "2020-03-25T21:05:39Z")
Thanks, but this was very buggy actually 😉 Here is a cleaned-up version and some results:
```julia
qq(xᵢ₋₂, xᵢ₋₁, xᵢ) = (xᵢ - xᵢ₋₁)/(xᵢ₋₁ - xᵢ₋₂)
function muller(P, oldest, older, old, atol=eps(one(old)))
@assert old ≠ older ≠ oldest ≠ old # we want q to be non-degenerate
xᵢ₋₂, xᵢ₋₁, xᵢ = oldest, older, old
# @show abs(xᵢ₋₂ - xᵢ)
while abs(xᵢ₋₁ - xᵢ) > atol
q = qq(xᵢ₋₂, xᵢ₋₁, xᵢ)
q² = q^2
q1 = q+one(q)
pᵢ = P(xᵢ)
pᵢ₋₁= P(xᵢ₋₁)
pᵢ₋₂= P(xᵢ₋₂)
A = q*pᵢ - q*q1*pᵢ₋₁ + q²*pᵢ₋₂
B = (q1+q)*pᵢ - q1^2*pᵢ₋₁ + q²*pᵢ₋₂
C = q1*pᵢ
den = let
Δ = sqrt(B^2 - 4A*C)
d⁺ = B + Δ
d⁻ = B - Δ
abs(d⁺) > abs(d⁻) ? d⁺ : d⁻
end
x = xᵢ - (xᵢ - xᵢ₋₁)*(2C/den)
xᵢ₋₂, xᵢ₋₁, xᵢ = xᵢ₋₁, xᵢ, x
# @show xᵢ₋₂, xᵢ₋₁, xᵢ, abs(xᵢ₋₁- xᵢ)
end
return xᵢ
end
f(x) = x^2-2
g(x) = x^3-1
using BenchmarkTools
@btime muller($f, 1.5, 1.0, 2.0) # → 1.41….; 42.341 ns (0 allocations: 0 bytes)
@btime muller($g, 1.5, 1.0, 2.0) # → 1.00; 231.045 ns (0 allocations: 0 bytes)
@time muller(f, 1.5, 1.0, 2.0) # → 1.41….
@time muller(f, 1.5, 1.0, -2.0) # → -1.41….
@time muller(g, 0.5, 0.5im, -0.5) # → -0.500 + 0.866…im
@time muller(g, 0.5, -0.5im, -0.5) # → -0.500 - 0.866…im
@time muller(g, -0.5, 0.0, 0.5) # → 1.00
```
unfortunately determinant is often negative (for random, real inputs) and throws `DomainError`. Any idea how to deal with this?
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 25, 2020, 9:19pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/35 "2020-03-25T21:19:18Z")
@abulak Maybe you could PR this to e.g. the `Roots.jl` package. cc @j_verzani
---
### Author:  [@hongchengni](https://discourse.julialang.org/u/hongchengni)
#### Post date: [March 25, 2020, 9:59pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/36 "2020-03-25T21:59:18Z")
I come across the Matlab program “Global complex Roots and Poles Finding algorithm” written by P. Kowalczyk:
> **[GitHub - PioKow/GRPF: Global complex Roots and Poles Finding algorithm](https://github.com/PioKow/GRPF)**
>
> Global complex Roots and Poles Finding algorithm. Contribute to PioKow/GRPF development by creating an account on GitHub.
It uses Delaunay triangulation to adaptively refine triangular mesh to locate all the roots (as well as poles). It looks like a very good algorithm.
It would be great if a translation into Julia could be done.
---
### Author:  [@abulak](https://discourse.julialang.org/u/abulak)
#### Post date: [March 26, 2020, 3:44am UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/37 "2020-03-26T03:44:20Z")
[https://github.com/JuliaMath/Roots.jl/pull/176](https://github.com/JuliaMath/Roots.jl/pull/176)
I’d be more interested to have a complex root finder in IntervalRootFinders, actually 😉
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 26, 2020, 3:53am UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/38 "2020-03-26T03:53:43Z")
> [@abulak](#):
>
> I’d be more interested to have a complex root finder in IntervalRootFinders, actually 😉
Me too… Would you like to help?
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [March 26, 2020, 7:32am UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/39 "2020-03-26T07:32:35Z")
It has an MIT license so no problem there. You should just go ahead and translate it to Julia!
---
### Author:  [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)
#### Post date: [May 5, 2020, 6:26pm UTC](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989/40 "2020-05-05T18:26:26Z")
Just came across this - I’ve actually already coded up the GRPF algorithm, but it’s not registered:
> **[GitHub - fgasdia/RootsAndPoles.jl: Julia implementation of the global complex...](https://github.com/fgasdia/RootsAndPoles.jl)**
>
> Julia implementation of the global complex root and pole finding (GRPF) algorithm. - GitHub - fgasdia/RootsAndPoles.jl: Julia implementation of the global complex root and pole finding (GRPF) algor...
Although I wouldn’t mind an Interval method either.
GRPF.jl works fine, but could probably use some performance improvements.
[Previous page](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989.md?page=1)
[Next page](https://discourse.julialang.org/t/complex-root-finder-for-a-general-function-f-z/35989.md?page=3)