Error with a Bordered Problem for a nonlinear system of 4 variables

[Oriane-Foussadier]

Hello, I'm new in the field of bifurcation analysis but now it's been a month than I'm passing through the documentation and different examples in order to apply the different (wonderful) methods you've implemented for a 1D periodic non-linear system of 4 variables with a global (integral) scalar condition linking two of my variables. I've seen that this condition can be implemented through the functional defined by the function BorderedProblem (my code is inspired by test-bordered-problem.jl). But I'm not able to make it work, even after several different things, and my last problem seems to be because of the use of BorderedArray which (if I understood well) must take our first guess as first component and the Lagrange multiplier as the last component, and which leads to a MethodError about the type for the extractVector.

I searched different things but I'm not able to solve this problem and not even sure how to continue (when this error will be solved) after for the continuation methods and finding the different branches etc.

Could you tell me where I'm wrong and indicate me briefly if there are some particular points for the rest of the procedure (continuation method, bifurcation points, branch switching)?

Thank you in advance for your help !

using Revise, ForwardDiff, Parameters, Setfield, Plots, LinearAlgebra
using Random
using Roots
using NumericalIntegration
using UnPack
using DiffEqOperators
using BifurcationKit, SparseArrays
const BK = BifurcationKit

#PARAMETERS
# All quantities are expressed in non-dimensional units
# Space grid
L  = 10*π                           # System size 
Nx = 512                            # Number of grid nodes
Δx = L / Nx                         # Grid spacing
x  = [-L/2 + i*Δx for i = 0:Nx-1]   # Space vector
# Model parameters (saved in a dictionary)
Pars=(Dc = 1e-2, Da = 1e-1, A = 1, Kd = 1, Ω0 = 0.5, Ω = 6.05, Ωd = 10, Z = 15, B = 15)

#Used norm in Newton
norminf(f) = norm(f, Inf)

function Field!(f, f0, Pars)
    @unpack Dc, Da, A, Kd, Ω0, Ω, Ωd, Z, B = Pars

    C  = @view f0[1:Nx]            
    Na = @view f0[Nx+1:2Nx]
    Ni = @view f0[2Nx+1:3Nx]
    V  = @view f0[3Nx+1:4Nx]

    Dx = CenteredDifference(1 ,3 , Δx, Nx) * PeriodicBC(Float64)
    Dxx= CenteredDifference(2 ,3 , Δx, Nx) * PeriodicBC(Float64)

    f[1:Nx]   = Dc * (Dxx * C) - (Dx * V).* C - (Dx * C).* V + A * Na - Kd * C
    f[Nx+1:2Nx] = Da * (Dxx * Na) - (Dx * V).* Na - (Dx *Na).* V + Ω0 * (1 .+ Ω .* Na.^2).* Ni - Ωd * C .* Na
    f[2Nx+1:3Nx]= (Dxx * Ni) - (Dx * V).* Ni - (Dx * Ni).* V - Ω0* (1 .+ Ω .* Na.^2).* Ni + Ωd.* C.* Na    
    f[3Nx+1:4Nx]= (Dx * C).* 2Z .* C - (Dx * C) .* 3B .* C .^ 2 + Dxx * V
    return f
end

Field(f0,Pars)=Field!(similar(f0),f0,Pars)

function Jacobian(f0, Pars)
    JacoField(f0)=Field(f0, Pars)
    return ForwardDiff.jacobian(JacoField, f0)
end

# Homogeneous Steady State (used in initial conditions below)
V0 = 0.
# Solving EqnHSS = 0 for Na gives Na at HSS
Eqn0(Na0) = - (Pars.Ωd * Pars.A / Pars.Kd) * Na0^2 + Pars.Ω0 * (1 - Na0) * (1 + Pars.Ω * Na0^2)

Na0   = find_zero(Eqn0, (0, 1)) #Find the roots of Eqn0 in the intervall 0,1

# Ni and C at HSS are readily calculated from NaHSS
Ni0 = 1 - Na0
C0  = (Pars.A / Pars.Kd) * Na0

# Guess for the stationary solution
sol0 = vcat(C0.* ones(Nx), Na0.* ones(Nx), Ni0.* ones(Nx), V0.* ones(Nx))

# Conservative condition on Na and Ni (integral over space of Na+Ni == 512)
g(f, Pars)= sum(f[Nx+1:3Nx])-512
probBordered = BorderedProblem(Field, g, (@lens _.Ω)) 

#Newton for our Bordered Problem
optnewton = NewtonPar(verbose = false)
optnewbd  = (@set optnewton.linsolver = MatrixBLS(optnewton.linsolver))
z0    = BorderedArray(sol0, 3.)
newtonBordered(probBordered, z0, Pars, optnewbd)

[rveltz]

Sorry for that but the BorderedProblem is really not developed as mentioned in the docs. For the moment, you should define a new problem where you append the g functional to Field.

Also, you should not compute Dx, Dxx at every call of Field but rather store them in Pars.

[rveltz]

Did you succeed with this?

# you should put Nx in Pars so that it is not global
function FF(f, f0, Pars)
    Field!(f, f0, Pars)
    f[end] = sum(f[Nx+1:3Nx])-512
    f
end

and then prob_bif = BifurcationProblem(FF, ...)

[Oriane-Foussadier]

Hello, yes the function "works" but according to our theoretical previsions it doesn't match. I need to check some properties of the jacobian before going further. As we have a conserved quantity in our system, we must have a null jacobian... which is not the case, and if so, we will certainly have a convergence problem in the Newton method. We plan to maybe re-express our system to find the roots of a system derived from the Lagrange multipliers. So... still in progress.

Maybe a last question to be sure, is it not possible to pass a non squared (I mean a non n*m matrix, n != m) Jacobian to the Bifurcation Problem? I tried but maybe the problem was elsewhere...

[rveltz]

I see. You should try newton with DefaultLS() and report the error here. Also, for continuation, you should use MoorePenrose() which basically is a newton solver on overparametrised equation. So this should (perhaps) handle your edge case.

But then the bifurcation detection will be brittle.

Multipliers seem a better.

Also, this looks like a DAE no?

[rveltz]

I read HSS, so you may be doing Slow-fast analysis? can you not put back the epsilon?

[Oriane-Foussadier]

Hello, sorry for the delay.

In our case HSS stands for Homogeneous Steady State... I'm not sure about what you're referring to by "Slow-fast analysis" but I think it is not what we are doing.

For the moment we left aside the Lagrange multipliers because we managed to have a non-zero jacobian by suppressing downright an equation for a field at one spatial point and reexpressing it with the constraint (so not an actual dynamical equation) but which can make it work (and also I had to implement the jacobian manually, I don't know yet why but the function ForwardDiff.jacobian was not relevant in this case). However we have difficulties to trigger the continuation steps/the tangent predictor because of either Newton algo that failed to converge or even more : a singular exception. So we are currently testing different values of parameters...

And yes we have a DAE and we have 3 fields with temporal dependence but also the velocity field with no temporal dependence. Do you have a special recommendation about that?

[rveltz]

I don't know yet why but the function ForwardDiff.jacobian was not relevant in this case).

PreallocationTools.jl might help there, but actually it seems it is not needed:

This works for me:

@views function Fvf!(f,f0, Pars)
    Field!(f[1:end-1], f0, Pars)
    f[end] = g(f0, Pars)
    f
end
Fvf(f,Pars)=Fvf!(similar(f,length(f)+1),f,Pars)

prob = BifurcationProblem(Fvf,sol0,Pars)

newton(prob, NewtonPar(verbose = true))

[Oriane-Foussadier]

Oh no sorry I misspoke I said that Newton failed to converged but I meant after only very a few continuation iterations, even with very small ds, and different values for theta. I put back our code here maybe you can see it directly (we've changed a little bit the representation of our equation (instead of real concentration Ni and Na, we have now Np=Ni+Na and Nm= Na-Ni, but it's basically the same):

function Field!(f, f0, Pars)
    @unpack Dc, Da_p, Da_m, A, Kd, Ω0, Ω, Ωd, Z, B, Dx, Dxx = Pars
    C  = @view f0[1:Nx]
    Np = @view f0[Nx+1:2Nx]
    Nm = @view f0[2Nx+1:3Nx]
    V  = @view f0[3Nx+1:4Nx]

    f[1:Nx] = Dc * (Dxx * C) - (Dx * V).* C - (Dx * C).* V + A * ( Np .+ Nm ) / 2 - Kd * C
    f[Nx+1] = Nx - sum(Np)
    f[Nx+2:2Nx] = 0.5 * (Da_p * (Dxx * Np)[2:end] .+ (Da_m * (Dxx * Nm)[2:end])) - ((Dx * V).* Np)[2:end] - ((Dx *Np).* V )[2:end]
    f[2Nx+1:3Nx]= 0.5 * (Da_m * (Dxx * Np) + (Da_p * (Dxx * Nm))) - (Dx * V).* Nm - (Dx *Nm).* V + Ω0* (1 .+ Ω .* (0.5*(Np + Nm)).^2).* (Np - Nm) - Ωd.* C.* (Np + Nm)   
    f[3Nx+1:4Nx]= (Dx * C).* 2Z .* C - (Dx * C) .* 3B .* C .^ 2 + Dxx * V
    return f
end
Field(f0,Pars)=Field!(similar(f0),f0,Pars)

# Homogeneous Steady State (used in initial conditions below)
V0  = 0.
Np0 = 1.
# Solving EqnHSS = 0 for Na gives Na at HSS
Eqn0(Na0) = - (Pars.Ωd * Pars.A / Pars.Kd) * Na0^2 + Pars.Ω0 * (1 - Na0) * (1 + Pars.Ω * Na0^2)
Na0   = find_zero(Eqn0, (0, 1)) #Find the roots of Eqn0 in the intervall 0,1
# Ni and C at HSS are readily calculated from NaHSS
Nm0 = 2*Na0-Np0
C0  = (Pars.A / Pars.Kd) * Na0

# Guess for the stationary solution
sol0 = vcat(C0.* ones(Nx), Np0.* ones(Nx), Nm0.* ones(Nx), V0.* ones(Nx))

probBif = BK.BifurcationProblem(Field, sol0, Pars, (@lens _.Ω) ;
  J = dField ,
  plotSolution = (f, Pars; kwargs...) -> (plotsol(f; label="", kwargs... )),
  recordFromSolution = (f, Pars; kwargs...) -> (Cmean = integrate(x,f[1:Nx])/L , field = f)
    )

opt_newton = NewtonPar(eigsolver=EigArpack(1.01, :LM), verbose=true, tol = 0.00001)

opts_br_eq = ContinuationPar(dsmin = 0.000000001, dsmax = 0.001, ds = 0.000001 , θ=0.5,
    pMin = 6.05 , pMax = 10. , detectBifurcation = 3, nev = 21,
    newtonOptions = opt_newton, maxSteps = 50000,
    nInversion = 20)

br = continuation(probBif, PALC(), opts_br_eq, verbosity=5, normC = norminf)

[rveltz]

OK, I am missing a few things, this is not runnable.

[Oriane-Foussadier]

Oh yes sorry the parameter are the same as before and the jacobian is implemented in another file as it is long. Here you have :

norminf(f) = norm(f, Inf)

L  = 10*π                           # System size   
Nx = 512                       # Number of nodes
Δx = L / Nx                       # Grid spacing
x  = [-L/2 + i*Δx for i = 0:Nx-1]   # Space vector
Pars=(Dc = 1e-2, Da_p = 1e-1+1, Da_m = 1e-1-1, A = 1, Kd = 1, Ω0 = 0.5, Ω = 6.05, Ωd = 10, Z = 15, B = 15, 
    Dx = CenteredDifference(1 ,3 , Δx, Nx) * PeriodicBC(Float64), 
    Dxx = CenteredDifference(2 ,3 , Δx, Nx) * PeriodicBC(Float64))

function dField!(f, f0, Pars)
     @unpack Dc, Da_p, Da_m, A, Kd, Ω0, Ω, Ωd, Z, B, Dx, Dxx = Pars
    C  = @view f0[1:Nx]            
    Np = @view f0[Nx+1:2Nx]
    Nm = @view f0[2Nx+1:3Nx]
    V  = @view f0[3Nx+1:4Nx]
#############################################################      Jc     ###############################################################################
   #Jcc
    diag_cc = ( -Dc / (2 * Δx ^2 )  .- Dx*V .- Kd )   .*   ones(eltype(f), Nx)
    m1_cc = ( V[2:end] / (2 * Δx) )                   .* ones(eltype(f), Nx-1)  
    m2_cc = ( Dc / (4 * Δx ^2 ) )                     .* ones(eltype(f), Nx-2) 
    p1_cc = -( V[1:end-1] / (2 * Δx) )                .* ones(eltype(f), Nx-1) 
    mn_cc = ( V[1] / (2 * Δx) )                      .* ones(eltype(f), 1)
    mn1_cc = (Dc / (4 * Δx^2) )                      .* ones(eltype(f), 2) 
    pn_cc = ( -V[end] / (2 * Δx) )                   .* ones(eltype(f), 1)

    Jcc=spdiagm(0 => diag_cc, 1 => p1_cc, 2 => m2_cc, -1 => m1_cc, -2 => m2_cc, -Nx+1 => pn_cc, -Nx+2 => mn1_cc, Nx-1 => mn_cc, Nx-2 => mn1_cc)

    #JcNp = JcNm
    diag_cNp   = (A/2) .* ones(eltype(f), Nx) 

    JcNp=spdiagm(0 => diag_cNp)
    JcNm=JcNp

    #JcV
    diag_cV  = (- Dx * C )                   .* ones(eltype(f), Nx)    
    m1_cV = ( C[2:end] / (2 * Δx) )          .* ones(eltype(f), Nx-1) 
    p1_cV = -( C[1:end-1] / (2 * Δx) )       .* ones(eltype(f), Nx-1)
    mn_cV = ( C[1] / (2 * Δx) )              .* ones(eltype(f), 1)
    pn_cV = -( C[end] / (2 * Δx) )           .* ones(eltype(f), 1)

    JcV=spdiagm(0 => diag_cV, -1 => m1_cV, 1=> p1_cV, Nx-1 => mn_cV, -Nx+1 =>pn_cV)

    ######JC
    JC=hcat(Jcc, JcNp, JcNm, JcV)
#############################################################      JNp    ###############################################################################
    #JNpc
    diag_Npc   = zeros(eltype(f), Nx) 

    JNpc=spdiagm(0 => diag_Npc)

    #JNpNp
    diag_NpNp = (-Da_p / (4* Δx^2) .- (Dx * V))                    .* ones(eltype(f), Nx)
    m1_NpNp = ( V[2:end] / (2 * Δx) )                              .* ones(eltype(f), Nx-1)  
    m2_NpNp = ( Da_p / (8* Δx^2) )                                 .* ones(eltype(f), Nx-2) 
    p1_NpNp = ( -V[1:end-1] / (2 * Δx) )                           .* ones(eltype(f), Nx-1)
    mn_NpNp = -1                                                   .* ones(eltype(f), 1)
    mn1_NpNp = ( Da_p / (8* Δx^2) )                                .* ones(eltype(f), 2)
    pn_NpNp = ( -V[end] / (2 * Δx) )                               .* ones(eltype(f), 1)

    JNpNp=spdiagm(0 => diag_NpNp, 1 => p1_NpNp, 2 => m2_NpNp, -1 => m1_NpNp, -2 => m2_NpNp, -Nx+1 => pn_NpNp, -Nx+2 => mn1_NpNp, Nx-1 => mn_NpNp, Nx-2 => mn1_NpNp)
    JNpNp[1,:].=-1

    #JNpNm
    diag_NpNm = (- Da_m / (4* Δx^2) )                                .* ones(eltype(f), Nx)
    m2_NpNm =   ( Da_m / (8* Δx^2) )                                 .* ones(eltype(f), Nx-2) 
    mn1_NpNm = ( Da_m / (8* Δx^2) )                                  .* ones(eltype(f), 2)

    JNpNm=spdiagm(0 => diag_NpNm, 2 => m2_NpNm, -2 => m2_NpNm, -Nx+2 => mn1_NpNm, Nx-2 => mn1_NpNm)
    JNpNm[1,:].=0

    #JNpV
    diag_NpV  = (- Dx * Np )                   .* ones(eltype(f), Nx)    
    m1_NpV = ( Np[2:end] / (2 * Δx) )          .* ones(eltype(f), Nx-1) 
    p1_NpV = -( Np[1:end-1] / (2 * Δx) )       .* ones(eltype(f), Nx-1)
    pn_NpV = -( Np[end] / (2 * Δx) )           .* ones(eltype(f), 1)

    JNpV=spdiagm(0 => diag_NpV, -1 => m1_NpV, 1=> p1_NpV, -Nx+1 =>pn_NpV)
    JNpV[1,:].=0

    ######JNp
    JNp=hcat(JNpc, JNpNp, JNpNm, JNpV)   
#############################################################      JNm    ###############################################################################
    #JNmc
    diag_Nmc   = -( Ωd * (Np.- Nm) )     .* ones(eltype(f), Nx) 

    JNmc=spdiagm(0 => diag_Nmc)

    #JNmNp
    diag_NmNp = (-Da_m / (4* Δx^2) .+ Ω0.*( (Np .- Nm) .* (0.5 .* Ω .* (Np .+ Nm)) .+ (1 .+ (Ω/4)*(Np .+ Nm).^2)) .- Ωd .* C)    .* ones(eltype(f), Nx)
    m2_NmNp = ( Da_m / (8* Δx^2) )                                 .* ones(eltype(f), Nx-2) 
    mn1_NmNp = ( Da_m / (8* Δx^2) )                               .* ones(eltype(f), 2)

    JNmNp=spdiagm(0 => diag_NmNp, 2 => m2_NmNp, -2 => m2_NmNp, -Nx+2 => mn1_NmNp, Nx-2 => mn1_NmNp)

    #JNmNm
    diag_NmNm = (-Da_p / (4* Δx^2) .+ Ω0.*( (Np.-Nm) .* (0.5 * Ω * (Np .+ Nm)) .- (1 .+ (Ω/4).*(Np .+ Nm).^2)) .- Ωd.* C .- Dx*V )    .* ones(eltype(f), Nx)
    m1_NmNm = ( V[2:end] / (2 * Δx) )                             .* ones(eltype(f), Nx-1)  
    m2_NmNm = ( Da_p / (8* Δx^2) )                                .* ones(eltype(f), Nx-2) 
    p1_NmNm = ( -V[1:end-1] / (2 * Δx) )                          .* ones(eltype(f), Nx-1)
    mn_NmNm = ( V[1] / (2 * Δx) )                                 .* ones(eltype(f), 1) 
    mn1_NmNm = ( Da_p / (8* Δx^2) )                               .* ones(eltype(f), 2)
    pn_NmNm = ( -V[end] / (2 * Δx) )                              .* ones(eltype(f), 1)

    JNmNm=spdiagm(0 => diag_NmNm, 1 => p1_NmNm, 2 => m2_NmNm, -1 => m1_NmNm, -2 => m2_NmNm, -Nx+1 => pn_NmNm, -Nx+2 => mn1_NmNm, Nx-1 => mn_NmNm, Nx-2 => mn1_NmNm)

    #JNmV
    diag_NmV  = (- Dx * Nm )                  .* ones(eltype(f), Nx)    
    m1_NmV = ( Nm[2:end] / (2 * Δx) )         .* ones(eltype(f), Nx-1)
    mn_NmV = ( Nm[1] / (2 * Δx) )             .* ones(eltype(f), 1) 
    p1_NmV = -( Nm[1:end-1] / (2 * Δx) )      .* ones(eltype(f), Nx-1)
    pn_NmV = -( Nm[end] / (2 * Δx) )          .* ones(eltype(f), 1)

    JNmV=spdiagm(0 => diag_NmV, -1 => m1_NmV, 1=> p1_NmV, Nx-1 => mn_NmV, -Nx+1 =>pn_NmV)

    ######JNp
    JNm=hcat(JNmc, JNmNp, JNmNm, JNmV)
#############################################################      JV    ###############################################################################  
    #JVC
    diag_Vc= ((Dx*C) .* (2 * Z .- 6 * B * C))                       .* ones(eltype(f), Nx)
    m1_Vc = -(1/(2 * Δx)) * (2 * Z * C .- 6 * B * C.^2)[2:end]      .* ones(eltype(f), Nx-1)
    p1_Vc = (1/(2 * Δx)) * (2 * Z * C .- 6 * B * C.^2)[1:end-1]     .* ones(eltype(f), Nx-1) 
    mn_Vc = -(1/(2 * Δx)) * (2 * Z * C[1] .- 6 * B * C[1] .^2)      .* ones(eltype(f), 1)
    pn_Vc = (1/(2 * Δx)) * (2 * Z * C[end] .- 6 * B * C[end] .^2)   .* ones(eltype(f), 1)

    JVC= spdiagm(0 => diag_Vc, -1 => m1_Vc, 1=> p1_Vc, Nx-1 => mn_Vc, -Nx+1 =>pn_Vc)

    #JVNp = JVNm
    diag_VN = zeros(eltype(f), Nx)
    JVN = spdiagm(0 => diag_VN)

    #JVV
    diag_VV = -( 1 / (2 * Δx ^2) )          .* ones(eltype(f), Nx)
    m2_VV = +( 1 / (4 * Δx ^2) )           .* ones(eltype(f), Nx-2)
    mn_VV = +(1/(4 * Δx^2))                .* ones(eltype(f), 1)

    JVV= spdiagm(0 => diag_VV, -2 => m2_VV, 2 => m2_VV, Nx-1 => mn_VV, -Nx+1 => mn_VV)

    ###JV
    JV=hcat(JVC, JVN, JVN, JVV)

    jacobian = vcat(JC, JNp, JNm, JV)
    return jacobian
end 

dField(f0, Pars) = dField!(similar(f0), f0, Pars)

[rveltz]

Ok, I found your mistake: you did not update the vector f in Ffield!!! note the .= Then AD works out of the box, although slowly because it produces a dense jacobian. you could improve it with sparity detection

function Field!(f, f0, Pars)
    @unpack Dc, Da_p, Da_m, A, Kd, Ω0, Ω, Ωd, Z, B, Dx, Dxx = Pars
    C  = @view f0[1:Nx]
    Np = @view f0[Nx+1:2Nx]
    Nm = @view f0[2Nx+1:3Nx]
    V  = @view f0[3Nx+1:4Nx]

    f[1:Nx] .= Dc .* (Dxx * C) .- (Dx * V) .* C .- (Dx * C) .* V .+ A .* ( Np .+ Nm ) ./ 2 .- Kd .* C
    f[Nx+1] = Nx - sum(Np)
    f[Nx+2:2Nx] .= 0.5 * (Da_p * (Dxx * Np)[2:end] .+ (Da_m * (Dxx * Nm)[2:end])) - ((Dx * V).* Np)[2:end] - ((Dx *Np).* V )[2:end]
    f[2Nx+1:3Nx] .= 0.5 * (Da_m * (Dxx * Np) + (Da_p * (Dxx * Nm))) - (Dx * V).* Nm - (Dx *Nm).* V + Ω0* (1 .+ Ω .* (0.5*(Np + Nm)).^2).* (Np - Nm) - Ωd.* C.* (Np + Nm)
    f[3Nx+1:4Nx] .= (Dx * C).* 2Z .* C - (Dx * C) .* 3B .* C .^ 2 + Dxx * V
    return f
end
Field(f0,Pars)=Field!(similar(f0),f0,Pars)

# Homogeneous Steady State (used in initial conditions below)
V0  = 0.
Np0 = 1.
# Solving EqnHSS = 0 for Na gives Na at HSS
Eqn0(Na0) = - (Pars.Ωd * Pars.A / Pars.Kd) * Na0^2 + Pars.Ω0 * (1 - Na0) * (1 + Pars.Ω * Na0^2)
Na0   = find_zero(Eqn0, (0, 1)) #Find the roots of Eqn0 in the intervall 0,1
# Ni and C at HSS are readily calculated from NaHSS
Nm0 = 2*Na0-Np0
C0  = (Pars.A / Pars.Kd) * Na0

# Guess for the stationary solution
sol0 = vcat(C0.* ones(Nx), Np0.* ones(Nx), Nm0.* ones(Nx), V0.* ones(Nx))

probBif = BK.BifurcationProblem(Field, sol0, Pars, (@lens _.Ω) ;
  # J = dField ,
  plotSolution = (f, Pars; kwargs...) -> (plot!(f; label="", kwargs... )),
  # recordFromSolution = (f, Pars; kwargs...) -> (Cmean = integrate(x,f[1:Nx])/L , field = f)
    )

opt_newton = NewtonPar(eigsolver=EigArpack(1.01, :LM), verbose=true, tol = 1e-9, maxIter = 10)

opts_br_eq = ContinuationPar(dsmin = 0.000000001, dsmax = 0.005, ds = 1e-3 , θ=0.5,
    pMin = 6.0 , pMax = 10. , detectBifurcation = 0, nev = 21,
    newtonOptions = opt_newton, maxSteps = 50000,
    nInversion = 20)

# you need to use the Secant predictor because the problem is not square
br = continuation(probBif, PALC(tangent = Secant()), opts_br_eq, verbosity=3, plot = true, normC = norminf)

[Oriane-Foussadier]

Thank you very much ! Indeed with the right adjustment of parameters it works perfectly well ! I've only two last questions. First, it's optional but when I activate the plotSolution I would like to have a graph with the 4 fields ploted in the same range (the real space range), not successively, with the scale of the last one on the left. For this I tried with a new helper plot function, but I obtain each time an error about subplots or others..., maybe you can see directly what I'm doing wrong :

function GRAPH!(field)
    C  = @view field[1:Nx]            
    Np = @view field[Nx+1:2Nx]
    Nm = @view field[2Nx+1:3Nx]
    V  = @view field[3Nx+1:4Nx]
    graph=Plots.plot(x, [C, Np, Nm], xlabel="Spatial repartition", labels=["C" "Np" "Nm"], right_margin = 17Plots.mm, legend=:topleft)
    graph=plot!(twinx(),x, V, label="V", linecolor="red", legend=:topright, formatter=:plain, right_margin = 17Plots.mm)
    return graph
end 

and then in ContinuationPar :

plotSolution = (field, Pars; kwargs...) -> (GRAPH!(field))

And my last and important question is about the saving of files that contains the informations of the branches, solutions etc... I saw the saveToFile=true option of the ContinuationPar , which create indeed a file (but I'm not really able to open it correctly) and in the same time gives me an error ... Is it actually possible to save every step (with the field, parameters, eigenvalues, an eventually with the plots shown during the continuation process) ?

[rveltz]

Great!

For the plot, you need to pass the keywoord argument.

function GRAPH!(field; kp...)
    C  = @view field[1:Nx]            
    Np = @view field[Nx+1:2Nx]
    Nm = @view field[2Nx+1:3Nx]
    V  = @view field[3Nx+1:4Nx]
    graph=Plots.plot!(x, [C, Np, Nm], xlabel="Spatial repartition", labels=["C" "Np" "Nm"], right_margin = 17Plots.mm, legend=:topleft, kp...)
    graph=plot!(twinx(),x, V, label="V", linecolor="red", legend=:topright, formatter=:plain, right_margin = 17Plots.mm, kp...)

    # you can also do:
    plot!(rand(10), subplot = 1)
end 
plotSolution = (field, Pars; kwargs...) -> (GRAPH!(field; kwargs...))

saveToFile saves the branch during computation. It saves 2 files. You are right, it does not save everything but you can improve it, it is thius function

function saveToFile(iter::AbstractContinuationIterable, sol, p, i::Int64, br::ContResult)
                if iter.contParams.saveToFile == false; return nothing; end
                filename = iter.filename
                # this allows to save two branches forward/backward in case
                # bothside = true is passed to continuation
                fd = iter.contParams.ds >=0 ? "fw" : "bw"

                # create a group in the JLD format
                jldopen(filename*".jld2", "a+") do file
                    if haskey(file, "sol-$fd-$i")
                        delete!(file, "sol-$fd-$i")
                    end
                    mygroup = JLD2.Group(file, "sol-$fd-$i")
                    mygroup["sol"] = sol
                    mygroup["param"] = p
                end

                jldopen(filename*"-branch.jld2", "a+") do file
                    if haskey(file, "branch"*fd)
                        delete!(file, "branch"*fd)
                    end
                    file["branch"*fd] = br
                end
            end

Is it actually possible to save every step (with the field, parameters, eigenvalues, an eventually with the plots shown during the continuation process) ?

The plot no! But if you pass your own finaliseSolution, you can do whatever you like (saving, ...)

[rveltz]

Should we close this?

[Oriane-Foussadier]

Yes sorry for the delay.