bifurcationkit / BifurcationKit.jl

A Julia package to perform Bifurcation Analysis
https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable
Other
303 stars 37 forks source link

Question: detect fold and branch point separately #160

Closed rdprins closed 3 months ago

rdprins commented 3 months ago

When generating a bifurcation diagram (1 param), both folds and branch points are give the type 'bp'. Is there a way to differentiate between them?

For more context, this is my finalise_solution function, which saves the results. I would like the variable TY to be either 'fold' or 'branch point' when applicable.

function append_file(filename, s)
    open(filename, "a") do file
        println(file,s)
    end
end;

function finalise_solution(z, tau, step, br, filename)
    if (length(br.specialpoint)>0) && (br.specialpoint[end].step == step)
        TY = br.specialpoint[end].type
        special_type=true
    else
        TY = "RG"
        special_type=false
    end

    if special_type || (step % 10 == 0)
        stability = br.stable[end]
        append_file(filename, join([step, z.p, TY, stability, z.u...], ";"))
    end

    return true
end;
rveltz commented 3 months ago

Hi,

Sorry, this is tricky at the moment as I dont use events to locate bifurcation (as MatCont do). If you can call get_normal_form(br, ind_fold), that would help you for example doing get_normal_form(br, ind_fold) isa BK.Fold

rdprins commented 3 months ago

Thank you for your quick reply!

Unfortunately it isn't possible to call get_normal_form(br, length(br.specialpoint), nev=nev) within finalise_solution, because an error occurs in NormalForms.jl Specifically, at line 80: λ = real(br.eig[bifpt.idx].eigenvals[bifpt.ind_ev]) indexing goes out of bounds.

BoundsError: attempt to access 45-element Vector{@NamedTuple{eigenvals::Vector{ComplexF64}, eigenvecs::Matrix{ComplexF64}, converged::Bool, step::Int64}} at index [46]

Would you know how to resolve this?

rveltz commented 3 months ago

Here is another way based on events. The idea is to locate Folds (using FoldDetectEvent) and bifurcations (using BifDetectEvent)

You will see points like foldbp because the event is a Fold and a bifurcation point at the same time. This could be resolved at a later stage by changing :foldbp to fold

I hope this is enough for you, sorry for the inconvenience. Because I want to be able to address bifurcations of kernel dimension >1 like in PDE, I do not use the event machinery of MatCont/Auto but just check for change in the number of unstable eigenvalues.

opts_br = ContinuationPar(p_min = -10.0, p_max = -0., dsmax = 0.1, n_inversion = 6, nev = 3, detect_event = 2, detect_bifurcation = 0, max_steps = 130, detect_fold = false)
br = continuation(prob, PALC(tangent = Bordered()), opts_br; 
        normC = norminf, 
        # event = BK.FoldDetectEvent,
        # event = BK.BifDetectEvent,
        event = BK.PairOfEvents(BK.FoldDetectEvent, BK.BifDetectEvent),
        bothside = false)
rveltz commented 3 months ago

Unfortunately it isn't possible to call get_normal_form(br, length(br.specialpoint), nev=nev) within finalise_solution, because an error occurs in NormalForms.jl Specifically, at line 80: λ = real(br.eig[bifpt.idx].eigenvals[bifpt.ind_ev]) indexing goes out of bounds.

Just to be sure. Can you call get_normal_form(br, ind) once the branch has been computed without your finalizer?

rdprins commented 3 months ago

Here is another way based on events. The idea is to locate Folds (using FoldDetectEvent) and bifurcations (using BifDetectEvent)

You will see points like foldbp because the event is a Fold and a bifurcation point at the same time. This could be resolved at a later stage by changing :foldbp to fold

I hope this is enough for you, sorry for the inconvenience. Because I want to be able to address bifurcations of kernel dimension >1 like in PDE, I do not use the event machinery of MatCont/Auto but just check for change in the number of unstable eigenvalues.

opts_br = ContinuationPar(p_min = -10.0, p_max = -0., dsmax = 0.1, n_inversion = 6, nev = 3, detect_event = 2, detect_bifurcation = 0, max_steps = 130, detect_fold = false)
br = continuation(prob, PALC(tangent = Bordered()), opts_br; 
        normC = norminf, 
        # event = BK.FoldDetectEvent,
        # event = BK.BifDetectEvent,
        event = BK.PairOfEvents(BK.FoldDetectEvent, BK.BifDetectEvent),
        bothside = false)

This works, thanks so much!

I have two additional questions:

In other words, can I always expect the same behaviour as when using bifurcationdiagram(PALC(tangent = Secant()), detect_bifurcation = 3, detect_event=0, detect_fold=true)?

rveltz commented 3 months ago

Do I understand correctly that BK.Bordered will at least be equally correct as BK.Secant, but might take some additional time to compute?

Yes it is a more precise tangent predictor, takes more time to compute than Secant. FoldDetectEvent does not work well with Secant predictor

rveltz commented 3 months ago

Will bifurcationdiagram() work exactly the same after setting detect_bifurcation = 0? (Which is needed because detect_bifurcation = 3 and detect_event = 2 aren't compatible). E.g., do we still use the bisection algorithm to calculate branch points (using e.g. the set values for n_inversion and max_bisection_steps)?

Yes, you can pass events to bifurcationdiagram the same way you do for continuation

rveltz commented 3 months ago

Ill close this issue. Re-open if you have any other problem.