5

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)
phipsgabler
  • 20,535
  • 4
  • 40
  • 60
tugberk21
  • 53
  • 4
  • You should try and post Minimal (non-)Working Examples that we can reproduce on our own systems. For example, your code (as it currently stands) contains a few errors that are totally unrelated to your question: `g` is never called, `a` and `Muller_method` are not defined, the `deflation` function definition lacks an `end` keyword... Having an reproducible example would help us help you! – François Févotte Sep 14 '20 at 19:36
  • Incidentally, it is often not a good idea to initialize untyped arrays with `[]`. If you know that `roots` will only ever contain floating-point numbers, it would be much better to initialize it like this, for example: `roots = Float64[]`. – François Févotte Sep 14 '20 at 19:39
  • I didn't share that part because I thought it wasn't wrong. I'm sorry for this was my first post. Thank you very much. – tugberk21 Sep 14 '20 at 20:38

1 Answers1

4

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
François Févotte
  • 19,520
  • 4
  • 51
  • 74