@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
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.
1 Like
Thanks for your help!
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:
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.
1 Like
looking at http://www.grad.hr/nastava/gs/prg/NumericalRecipesinC.pdf (p. 371) the Muller method seems to be
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?
@abulak Did you translate their code? Numerical recipes has a terribly restrictive software license and it’s basically not allowed to use their code.
1 Like
@dpsanders there is no code for this method in the book quoted; this is literal translation of formulas
OK great, thanks!
Very nice code, by the way 
@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
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.)
1 Like
sorry, qq got pasted on the 0-th line; i corrected now
Thanks, but this was very buggy actually Here is a cleaned-up version and some results:
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?
1 Like
I come across the Matlab program “Global complex Roots and Poles Finding algorithm” written by P. Kowalczyk:
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.
2 Likes
https://github.com/JuliaMath/Roots.jl/pull/176
I’d be more interested to have a complex root finder in IntervalRootFinders, actually
1 Like
Me too… Would you like to help?
It has an MIT license so no problem there. You should just go ahead and translate it to Julia!
Just came across this - I’ve actually already coded up the GRPF algorithm, but it’s not registered:
Although I wouldn’t mind an Interval method either.
GRPF.jl works fine, but could probably use some performance improvements.
5 Likes