Pitchfork branching code weirdness

[antoine-levitt]

MWE:

using PseudoArcLengthContinuation, LinearAlgebra, Plots
const PALC = PseudoArcLengthContinuation
using ForwardDiff
using Parameters, Setfield

function f(x, p)
    [p.α*x[1] + x[1]^3 ]
end
function Jf(x, p)
    Jf = zeros(1, 1)
    Jf[1] = p.α+3x[1]^2
    dx -> Jf*dx
end
Fmit = f

# eigensolver
struct MyStruct <: PALC.AbstractEigenSolver
end
function (m::MyStruct)(J, nev)
    return [(J(1))[1]], [1.], true, 1
end

function D(f, x, p, dx)
    ε=1e-7
    (f(x + ε*dx, p) - f(x - ε*dx, p))/(2ε)
end
d1Fmit(x,p,dx1) = D((z, p0) -> Fmit(z, p0), x, p, dx1)
d2Fmit(x,p,dx1,dx2) = D((z, p0) -> d1Fmit(z, p0, dx1), x, p, dx2)
d3Fmit(x,p,dx1,dx2,dx3) = D((z, p0) -> d2Fmit(z, p0, dx1, dx2), x, p, dx3)
jet = (f, d1Fmit, d2Fmit, d3Fmit)
jet = (f, Jf, d2Fmit, d3Fmit)

alpha_min = -10.
alpha_max = 10.
alpha_start = 2.

# options for Newton solver
opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 20, eigsolver = MyStruct())
opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 20, eigsolver = MyStruct(), linsolver=GMRESKrylovKit())

# options for continuation
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = -0.01, pMax = alpha_max, pMin = alpha_min,
	                  detectBifurcation = 2, nev = 30, plotEveryNsteps = 10, newtonOptions = opt_newton,
	                  maxSteps = 100, precisionStability = 1e-6, nInversion = 4, dsminBisection = 1e-7, maxBisectionSteps = 25)

sol_start, _, _ = newton( f, [1.], (α=.1,), opt_newton)
par_mit = (α=alpha_start,)

br, _ = @time PALC.continuation(
    Fmit, Jf,
    sol_start, par_mit, (@lens _.α), opts_br;
    # printSolution = (x, p) -> norm(x),
    # plotSolution = (x, p; kwargs...) -> plotsol!(x ; kwargs...),
    plot = false, verbosity = 3)

br2, _ = continuation(jet..., br, 1, setproperties(opts_br, ds = -0.01, maxSteps = 20), plot = false, printSolution = (x,p) -> x[1], verbosity=3, Jt=Jf)
br3, _ = continuation(jet..., br, 1, setproperties(opts_br, ds = 0.01, maxSteps = 20), plot = false, printSolution = (x,p) -> x[1], verbosity=3, Jt=Jf)

Plots.plot([br,br2,br3], marker=2)

Also plot just [br, br2] and [br, br3].

The pitchfork is well detected, and the two branches get explored. However,

1. the base of the two branches is at the same point, whereas I would have expected both branches to have opposite starting points

2. the starting point is off the horizontal axis by an amount relatively large compared to my ds, whereas I would have expected the starting point to be exactly on the horizontal axis. I guess this is easy to fix by just adding the bifurcation point to the branch.

Are both these expected?

[rveltz]

1. This is expected because both options give the same starting point on the bifurcated branch. To find this first point, the sign of ds is irrelevant because a pitchfork is detected. Then, how you travel on the bifurcated branch is determined by ds. If you want to find the branch with opposite sign, you could use the argument amp_factor

2. I am not sure I understand. I have

6-element Array{Float64,1}:
 -0.00978636554320561
  0.09892606540057562
  0.0                
  0.01               
  0.5                
  0.0  

-0.0097 seems close to 0.01 to me. You are right that ds does not tell you the parameter shift for the first point though, it is (a bit) more complicated than that. You can force the parameter value of the first point on the branch using the argument δp though

[rveltz]

For 2. you mean vertically? In which case the solution is sqrt(ds) = 0.1 which is "large". I could perhaps overwrite the ds to stay close to zero. Note that evaluating the size of the neighborhood in which the normal form is "valid" is possible but...

[antoine-levitt]

1 It still feels strange. I would expect the starting point to be on the horizontal axis, and only after a ds step do we travel on the two branches.

2 But it should be ds vertically, right? The curve is vertical there, so the arclength is the same as the vertical. Maybe I'm confused by something, but I would expect (barring ds adaptation) to see something like a parabola equally spaced. Instead the starting point is much too far away. I think the problem is in pnew = bp.p + abs(ds) * dsfactor, which is going to result in a step sqrt(ds) away in arclength. It should be something like ds^2, no?

[rveltz]

1. You want the first point of the bifurcated branch to be the bifurcation point (in your case (0,0))? I agree it would be better. This is a no brainer to add

2. The first step is pure Newton actually. So pnew = bp.p + abs(ds) * dsfactor = ds = -0.01 in your case. Then, it tries to solve f=0 which gives x=sqrt(0.01). So the first point on the bifurcated branch is (-abs(ds), sqrt(abs(ds))). Then continuation` kicks in

[antoine-levitt]

2 But then that's not a step of ds in arclength, it's a much bigger step, no?

[rveltz]

2 Well it is "pseudo" arc length but that just evades the question :D ... There is no point on the branch yet (except the bifurcation point which I discarded up to now). Are you confused by what the docs say or something else? In the method continuation, the "same" is done for the first two points, see here

[rveltz]

Ok, I got it. You expect the first point to be at distance ds of the bifurcation point and not sqrt(ds). That's fair and I will change this.