Only trivial solution in BifurcationKit

Specific Domains Numerics
, ,
1

I try to investigate pitchfork bifurcation, using the BifurcationKit. My problem looks as follow. I start from some definitions:

ND = 500
dx = 1.0/ND
x = collect(1:ND)*dx.-dx/2

J = +1.0
beta0 = 0.1
u0 = 0.01*ones(Float64,ND)
W = ones(Float64,ND,ND);

function rhs(u, beta)
    global J,dx,W
    return u - tanh.(J*beta*dx*W*u)
end

Then, I introduce bifurcation problem:

prob = BifurcationKit.BifurcationProblem(rhs, u0, beta0, 1)

and use

opt_newton = BifurcationKit.NewtonPar(tol=1e-9)
opts = BifurcationKit.ContinuationPar(p_min=0.0, p_max=5.0, max_steps=500, dsmax=2.e-2, dsmin=1.e-4, ds=1.e-4, nev=ND, newton_options=opt_newton);

I run continuation method as follows:

br = BifurcationKit.continuation(prob, BifurcationKit.PALC(),opts, bothside = true, normC=norm, verbosity=0)

I see the correct bifurcation point at p\approx 1 and extract it:

p_test = br.specialpoint[2].param

Then, I try to switch branch:

test = BifurcationKit.continuation(br, 2, opts; δp = -0.1)

My final goal is to obtain 2 vectors p_values and solution_u. I have tried to extract this solution as br.sol or test.sol. But everywhere I see only the (seems that) trivial solution. For instance:

Безымянный
Безымянный1266×282 16.1 KB

I have tried to manipulate with record_from_solution, but it does not give the desired result. I have tried to add:

record_from_solution = (u,beta; k...) -> (n2=norm(u),solu=u)

I have read the documentation, but I still do not understand properly what did I wrong.

2

It seems too far from the bifurcation point. This works:

 test = BifurcationKit.continuation(br, 2, opts; δp = -0.01)

I am surprised it does not converged though.

3

If you use debug mode ENV["JULIA_DEBUG"] = BifurcationKit, your code gives:

┌ Debug: The guess for the amplitude of the first periodic orbit on the bifurcated branch obtained by the predictor is not small: 3.8685530832137087. This may lead to convergence failure of the first newton step or select a branch far from the Hopf point.
│ You can either decrease `ds` or `δp` (which is  how far from the bifurcation point you want the branch of periodic orbits to start). Alternatively, you can specify a multiplicative factor `ampfactor` to be applied to the predictor amplitude.
...

You have the indication You can either decrease ds or δp albeit the message is wrong (it is for Hopf bifurcation and I will correct this).

4

I definitely need to obtain more intution about \delta p… My final goal is to obtain array of parameter values p_values and the mean of solution, i.e. mean(sol_vector), where sol_vector represents the solution for given parameter value p. What should I do for it?

My naive guess was that I can do it by single run of BifurcationKit.continuation, but now I am not sure

5

If the library were optimal, you should not have to worry too much about it…

6

You have several options:

1/ plot(test.param, [mean sol.x for sol in test.sol])

2/ pass record_from_solution = (u,beta; k...) -> (n2=norm(u),solu=u, mean = mean(x)) and the plot(test.param, test.mean)

7

But with record_from_solution I have to run first the continuation, detect the bifurcation points, then switch the branch. Is it correct?

8

yes, the record function is part of prob.

9

For my specific goals, the following scheme gives the desired result:

  1. With the given prob = BifurcationKit.BifurcationProblem(...) (which is similar to the provided problem up to specific details), run
    diagram = BifurcationKit.bifurcationdiagram(prob, BK.PALC(), 2, opts, bothside=true, normC=norminf,
        verbosity=0, plot=false, halfbranch=true)
  1. Having stored diagram object in memory, extract information about all possible branches:
    for i in eachindex(diagram.child)
        parameter_values = sort(diagram[i].γ.param)
        solution_norm = sort(diagram[i].γ.x)
    end

Using this scheme, I can run one function bifurcationdiagram one time, avoiding the second run for branch switching.

The field child (as I understand) has the same length as the number of detected bifurcations, so diagram[i].γ gives the data about i-branch and then diagram[i].γ.param represents the corresponding parameter values, while diagram[i].γ.x gives the solutions norm. I use sort, but it is not necessary. It seems that sometimes extracted information a little bit messy, but it is easy to handle with it.

Also, bifurcation points for i-th branch can be extracted as diagram[i].γ.bp.p. Finally, the solutions (not their norm) are available from diagram[i].γ.sol. It is the array of named tuples. Filtering procedures are useful to extract specific elements of this array, for instance:

filtered = filter(elem -> elem.p > bp, diagram[i].γ.sol)

will extract elements that corresponds to parameter values p larger than bp.

All the information about parameters and solution is stored in field γ of diagram object. In my specific example, manipulating this object was the simplest possible way to store bifurcating solutions.

I am not sure that my way is completely correct from the logic of BifurcationKit library, however I obtain the desired data. From my point of view, it will be useful to add some information about diagram = BifurcationKit.bifurcationdiagram(...) structure into the package documentation.

To investigate diagram structure, it is enough to use fieldsname(typeof(diagram)) and read the respective documentation about specific fields (for instance, specialpoints and so on).

10

It does not! diagram[i].γ[1] is a named typle containing, among other things, the return value of your record_from_solution .

Also, bifurcation points for i -th branch can be extracted as diagram[i].γ.bp

actually, it is:

from(diagram[i])

We probably need a getter like get_solution(br, i) for this one:

diagram[i].γ.sol

I am not sure that my way is completely correct from the logic of BifurcationKit library, however I obtain the desired data.

In general, I would say no because it is not robust. That you had to do this is a sign of bad API and you should open an issue. It is best to use getters because the field of internal structs can change.

From my point of view, it will be useful to add some information about diagram = BifurcationKit.bifurcationdiagram(...) structure into the package documentation.

You can draft a PR if you have an idea, even if it is not complete.

Thanks a lot for your return on the use case.

11

But what if I do not specify any record_from_solution? I have extracted solutions vector from diagram[i].γ.sol in for loop and have checked that the this gives exactly the same as diagram[i].γ.x. I understand your point with “getters”, actually I just haven’t noticed get_branch function. I will try to reproduce the same data with get_branch and will return.

12

If you do not specify record_from_solution, the BifurcationProblem uses record_sol_default which returns (x = norm(x),)

13

So, despite I have proposed the wrong schemde since I have not specified record_from_solution the code works as I expected. get_branch does not simplify data extraction: I still need to manipulate with diagram[i].γ. I have marked your response as the solution. I am sorry for my poor knowledge of documentation. Thank you for your fast responses!

14

There is nothing wrong, you just picked the relevant key (.x) I guess by inspecting br[1].

I still need to manipulate with diagram[i].γ

Yes, I need to provide a getter for this. You can open an issue on BifurcationKit if you want