Deflation Method for Multiple Root Finding

Question

I'm trying to implement the deflation method for finding multiple roots of a polynomial on Julia. In the deflation method, new functions are generated from the actual function divided by x - x_roots[i], and the root of new generated function must be found. But g(x) = g(x) / (x - x_root) gives me a Stack Overflow error, as probably it generated an infinite recursive relation. How may I generate a new function in each step?

function deflation(f::Function, order)
    roots=[]
    n=1
    g(x)=f(x)
    x_root=Muller_method(f,-5,0,5,1e-5,50)
    while n<=order
        x_root=Muller_method(a,-5,0,5,1e-5,50)
        g(x)=g(x)/(x-x_root)
        append!(roots,x_root)
        n=n+1
    end
return (roots)

Answer 1

Something like this causes infinite recursion:

julia> g(x) = x
g (generic function with 1 method)

julia> g(1)
1

julia> g(x) = g(x) / 2
g (generic function with 1 method)

julia> g(1)
ERROR: StackOverflowError:
Stacktrace:
 [1] g(::Int64) at ./REPL[3]:1 (repeats 79984 times)

This is because function (or method) definitions do not work like variable assignment: each re-definition of g(x) overwrites the previous one (note how g above only ever has one method). When a method definition refers to itself, it is recursion, i.e. it refers to its own version at the time when the function gets called.

As for your deflation method, perhaps you could define a new function that closes over the vector of currently found roots. Consider the following example to see how things would work:

julia> f(x) = (x-1) * (x-2)
f (generic function with 1 method)

julia> roots = Float64[]
Float64[]


# g is defined once, and accounts for all currently found roots
julia> g(x) = f(x) / reduce(*, x-root for root in roots; init=one(x))
g (generic function with 1 method)


# No roots are identified at the beginning
julia> g(1+eps())  
-2.2204460492503126e-16

julia> g(2+eps())
0.0


# But the results produced by g update as roots are identified
julia> push!(roots, 1.)
1-element Array{Float64,1}:
 1.0

julia> g(1+eps())
-0.9999999999999998

julia> g(2+eps())
0.0