recommendations on basic usage of package

[wsshin]

(Starting an issue as recommended by https://discourse.julialang.org/t/ann-new-version-of-package-bifurcationkit-jl-v0-1-11/75569/8)

I'm new to BifurcationKit.jl, or bifurcation analysis in general. Still, from the package description, I feel that the package could solve my problem, but I am having difficulty in making it work.

My understanding is that BifurcationKit.jl solves $F(u,p) = 0$ with a varying parameter $p$. Here, $F$ is a vector-valued function in general, and from the documentation I see that it is usually a discretized version of some differential equation.

My problem is much simpler than this general situation. My function $f(β, k)$ is scalar-valued. The solution $β$ and the parameter $k$ are also scalar. I already have $f(β, k)$, such that I can evaluate $f(β, k)$ for arbitrary $β$ and $k$. There could be multiple solutions $β$ for a given $k$, and the number of solutions can change with $k$. Suppose we have $n$ solutions $β_1(k_0), ..., β_n(k_0)$ for some $k_0$. My goal is to track the evolution of each of these solutions while varying $k$ from $k_0$.

Currently I am solving $f(β, k) = 0$ for each $k$ by Root.jl, and here is the result: The plotted numbers indicate the number of the solution $β$'s for each value $k$ of the horizontal axis. It is hard to distinguish multiple $β$'s in this plot, but if we plot $b$ and $V$, which are normalized versions of $β$ and $k$, they become distinguishable:

I attempted to generate the top curve out of the three in the above plot with BifurcationKit.jl. I followed this example, and the overall structure of the code looks as follows:

using BifurcationKit

f(β, p) = ...

kmin, kmax = ...
βmin, βmax = ...

p₀ = (..., k=kmax)
β₀ = maximum(find_zeros(β->f(β,p₀), βmin, βmax))

optnewton = NewtonPar(tol = 1e-11, verbose = true)
prob = BifurcationProblem((β,p)->[f(β[1],p)], [β₀], p₀, (@lens _.k),
                          plotSolution = (β, p; kwargs...) -> plot!(β; ylabel="solution", label="", kwargs...))
sol = @time newton(prob, optnewton)

and running the code generates

┌─────────────────────────────────────────────────────┐
│ Newton Iterations      f(x)      Linear Iterations  │
├─────────────┬──────────────────────┬────────────────┤
│       0     │       1.1726e-17     │        0       │
└─────────────┴──────────────────────┴────────────────┘
  0.590592 seconds (985.27 k allocations: 61.903 MiB, 4.87% gc time, 98.68% compilation time)
NonLinearSolution{Vector{Float64}, BifurcationProblem{BifFunction{var"#23#26", BifurcationKit.var"#6#21"{var"#23#26"}, Nothing, BifurcationKit.var"#4#19"{var"#23#26"}, Nothing, BifurcationKit.var"#9#25"{BifurcationKit.var"#d1Fad#23"{var"#23#26"}}, BifurcationKit.var"#11#27", BifurcationKit.var"#13#29", BifurcationKit.var"#15#31", Bool, Float64}, Vector{Float64}, NamedTuple{(:ν, :k), Tuple{Int64, Float64}}, Setfield.PropertyLens{:k}, var"#24#27", typeof(BifurcationKit.recordSolDefault)}, Vector{Float64}, Int64}([1.8440990245986763e7], ┌─ Bifurcation Problem with uType Vector{Float64}
├─ inplace:  false
├─ Symmetric: false
└─ Parameter: k, [1.1726324121788534e-17], true, 0, 0)

The solution was found in one Newton iteration, because I used the solution found by Roots.jl as the initial guess.

Then, if I continue and run

optcont = ContinuationPar(dsmin=0.01, dsmax=0.2, ds=-0.1, pMin=kmin, pMax=kmax,
	                      newtonOptions=NewtonPar(maxIter=10, tol=1e-9))

br = continuation(prob, PALC(), optcont; plot = true)

I get

 ┌─ Number of points: 101
 ├─ Curve of EquilibriumCont
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter k starts at 1.2566370614359174e7, ends at 1.2566354965157995e7
 ├─ Algo: PALC
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1, endpoint at k ≈ +12566354.80804161,                                                                     step = 101

Here, the endpoint k value is slightly moved from kmax = 12566370.614359174 towards kmin = 897597.9010256552 (because ds < 0), but it is much greater than kmin, whereas the smallest k value plotted in the top curve of the first figure above is actually quite close to kmin.

What is the right approach to solve this problem with BifurcationKit.jl? Any suggestions? If desired, I can share the full script used to generate the above result.

[rveltz]

I suggest that you first use Natural() continuation to get à feeling. I think the problem is from the scaling of k. I would write k = 1e6 * k2 and continue in k2

[wsshin]

Now my function accept k / 1e6 and convert it to k internally. Further, I tried Natural(), but it doesn't work:

┌─ Number of points: 2
 ├─ Curve of EquilibriumCont
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter k starts at 12.565703947692505, ends at 12.566370614359172
 ├─ Algo: Natural
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1, endpoint at k ≈ +12.56637061,                                                                     step =   1

If I use PALC(), then it generates the following erro:

ERROR: AmosException with id 2: overflow.
Stacktrace:
  [1] _besseli(nu::Float64, z::ComplexF64, kode::Int32)
    @ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:240
  [2] besseli(nu::Float64, z::ComplexF64)
    @ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:364
  [3] besseli
    @ ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:454 [inlined]
...

I think this is an error related to automatic differentiation. My function $f(β,k)$ is defined with Bessel functions in the complex domain; even though $β$ and $k$ are real, the arguments of the Bessel functions used in the definition of $f(β, k)$ become complex in some cases. This makes automatic differentiation by ForwardDiff.derivative() difficult.

Is there a way to pass the Jacobian function $∂f/∂β$ explicitly? (Note that my $f$ is scalar-valued, so $∂f/∂β$ is the full Jacobian.) I have a way to evaluate $∂f/∂β$ using the properties of the Bessel functions, without using ForwardDiff.derivative(). I would like to see if passing the externally constructed Jacobian suppresses the use of ForwardDiff.derivatine(), thereby removing the above error.

[rveltz]

You can pass the jacobian to a BifurcationProblem with ; J = ...)

[rveltz]

I am surprised Natural did not work, it is a basic Newton method

[rveltz]

Oh are using complex numbers? You can only use reals

[wsshin]

My function uses complex numbers internally, but its arguments $(β,k)$ and output $f$ are all real. This should be OK, right?

Providing the externally constructed Jacobian doesn't help unfortunately...

I increased the lower bound of $k$ from kmin = 2π/6 to kmin = 2π/3 to make the problem easier, and now I get the following result (with PALC; Natural still doesn't work)

 ┌─ Number of points: 67
 ├─ Curve of EquilibriumCont
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter k starts at 12.566370614359172, ends at 2.228729830345902
 ├─ Algo: PALC
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1, endpoint at k ≈ +2.09439510,                                                                     step =  66

The end point value k ≈ +2.09439510 is close to kmin = 2π/3 ≈ 2.0943951023931953, so I guess the algorithm is working, but I would like see the plot of $β(k)$ to verify if the solution is correct. How can I evaluate $β(k)$?

[rveltz]

You can look at br.branch or do plot(br) to check Beta(k).

I am surprised natural does not work. You should check the function newton in verbose mode. Once this works try Natural by adjusting the ds. Then try PALC. Do this with verbosity =3 in Newton verbose mode. This will tell you the parameters steps and the Newton iterations.

[wsshin]

Thanks. I retrieved $β$ and $k$ from br.branch and normalized them into $b$ and $V$. Here is the plot of $b$ vs. $V$: The plot replicates the corresponding plot in the original posting, so BifurcationKit.jl is working to some degree!

But note that the curve does not go all the way towards the origin (0,0). This is because I set the lower bound of k2 rather large, at k2min = 2π / 3.0, corresponding to kmin = k2min * 1e6. (By the way, the plot in the original posting was generated for k2min = 2π / 6.0.) If I decrease k2min even slightly, I get the following result:

Any suggestions for overcoming this difficulty? I tried to decrease the magnitudes of ds, 'dsmin', 'dsmax', but it didn't help much.

[wsshin]

Regarding Natural, here is the output with verbosity = 3, but I'm not sure how to modify the code based on this result:

#####################################################
────────── Natural ────────────

─────────────────  INITIAL GUESS ────────────────────
--> convergence of initial guess = OK

--> parameter = 12.566370614359172, initial step

───────────────── INITIAL TANGENT ───────────────────
--> convergence of the initial guess = OK

--> parameter = 12.565703947692505, initial step (bis)
──────────────────────────────────────────────────────────────────────
Continuation Step 0 
Step size = -1.0000e-01
Parameter k = 1.2566e+01 ⟶  1.2566e+01 [guess]
--> Step Converged in 0 Nonlinear Iteration(s)
Parameter k = 1.2566e+01 ⟶  1.2566e+01
--> Computed 1 eigenvalues in 1 iterations, #unstable = 0

In comparison, here is the abbreviated output for PALC:

#####################################################
────────── PALC ────────────

─────────────────  INITIAL GUESS ────────────────────
--> convergence of initial guess = OK

--> parameter = 12.566370614359172, initial step

───────────────── INITIAL TANGENT ───────────────────
--> convergence of the initial guess = OK

--> parameter = 12.565703947692505, initial step (bis)
Predictor:  Secant
──────────────────────────────────────────────────────────────────────
Continuation Step 0 
Step size = -1.0000e-01
Parameter k = 1.2566e+01 ⟶  1.2488e+01 [guess]
--> Step Converged in 1 Nonlinear Iteration(s)
Parameter k = 1.2566e+01 ⟶  1.2488e+01
--> Computed 1 eigenvalues in 1 iterations, #unstable = 0
Predictor:  Secant
──────────────────────────────────────────────────────────────────────
Continuation Step 1 
Step size = -1.4050e-01
Parameter k = 1.2488e+01 ⟶  1.2377e+01 [guess]
--> Step Converged in 1 Nonlinear Iteration(s)
Parameter k = 1.2488e+01 ⟶  1.2377e+01
--> Computed 1 eigenvalues in 1 iterations, #unstable = 0
Predictor:  Secant
──────────────────────────────────────────────────────────────────────
...
──────────────────────────────────────────────────────────────────────
Continuation Step 66 
Step size = -2.0000e-01
Parameter k = 2.2287e+00 ⟶  2.0874e+00 [guess]
Newton correction failed
Halving continuation step, ds=-0.1
Predictor:  Secant

The message for the last continuation step seems to indicate a failure. Any suggestions for preventing this failure?

[rveltz]

suggestions for overcoming this difficulty?

Yes I think you should decrease the option maxIter, somehow Newton works too well. Indeed it found a solution far from the guess as you can see with the jump

[rveltz]

Regarding Natural, here is the output with verbosity = 3, but I'm not sure how to modify the code based on this result:

I am wondering what happens if pmax is a bit higher than your initial guess? I don't get why the continuation ends. Can you show the Newton iterations too?

I am sorry about this, being away from my comp makes it hard to help you. Usually this would be resolved quickly

[wsshin]

Can you show the Newton iterations too?

Could you let me know how to show the Newton iterations? The posted output is all I got with verbosity = 3.

[rveltz]

That is in NewtonPar

[rveltz]

can you share the function F please so I can help more?

[rveltz]

can we progress on this?

[wsshin]

Sorry for the delay. I think I found a way to solve my equation without using this package. The method modified my function, and I think there is a high chance that the modified function could work better with BifurcationKit.jl as well.

Even though I have found a solution not using BifurcationKit.jl, I am curious if the modified function would work better with the package than the previous function. I will perform a test and post the result when I have a chance. Until then, I will close this issue. Thanks for your help so far!

[rveltz]

Please do! I havent tried examples with large values like you do.

Bests