Skip to content

Instantly share code, notes, and snippets.

@rveltz

rveltz/catalystexp1.jl Secret

Created May 12, 2021
Embed
What would you like to do?
using Revise
using BifurcationKit, Plots, LinearAlgebra, Setfield
const BK = BifurcationKit
using Catalyst
rn = @reaction_network begin
(v0+v*(S*X)^n/((S*X)^n+(D*A)^n+K^n),d), ∅ ↔ X
*X,τ), ∅ ↔ A
end S D τ v0 v K n d
odefun = ODEFunction(convert(ODESystem,rn),jac=true)
F = (u,p) -> odefun(Array(u),p,0) # this is because MK does not like views
J = (u,p) -> odefun.jac(Array(u),p,0) # this is because MK does not like views
params = [1.,9.,0.001,0.01,2.,20.,3,0.05]
p_idx = 1 # The index of the bifurcation parameter.
p_span = (0.1,20.) # The parameter range for the bifurcation diagram.
plot_var_idx = 1 # The index of the variable we plot in the bifurcation diagram.
opts = ContinuationPar( dsmax = 0.95, # Maximum arclength value of the pseudo-arc length continuation method.
dsmin = 1e-4, # Minimum arclength value of the pseudo-arc length continuation method.
ds=0.001, # Initial arclength value of the pseudo-arc length continuation method (should be positive).
maxSteps = 130, # The maximum number of steps.
pMin = p_span[1], # Minimum p-vale (if hit, the method stops).
pMax = p_span[2], # Maximum p-vale (if hit, the method stops).
detectBifurcation = 3, # Value in {0,1,2,3} determining to what extent bifurcation points are detected (0 means nothing is done, 3 both them and there localisation are detected).
nInversion = 6, nev = 2,
newtonOptions = NewtonPar(tol = 1e-9, verbose = false, maxIter = 5)) #Parameters to the newton solver (when finding fixed points) see BifurcationKit documentation.
params_input = setindex!(copy(params),p_span[1],p_idx) # The input parameter values have to start at the first index of our parameter span.
@error 1==0 "Ca ne marche pas du tout sur les valeurs propres!!!"
br, = continuation(F, J, fill(0.,length(rn.states)), params_input, (@lens _[p_idx]) ,opts; # Gives our input.
linearAlgo = MatrixBLS(),
# tangentAlgo = BorderedPred(),
verbosity = 3, plot=true, # We do not want to display, or plot, intermediary results.
printSolution=(x, p) -> x[plot_var_idx], # How we wish to print the output in the diagram. Here we simply want the value of the target varriable.
# perturbSolution = (x,p,id) -> (x .+ 0.8 .* rand(length(x))), # Parameter for the continuation method, see BifurcationKit documentation.
# callbackN = (x, f, J, res, iteration, itlinear, options; kwargs...) -> res <1e7 # Parameter for the continuation method, see BifurcationKit documentation.
)
####################################################################################################
# computation of the derivatives
using ForwardDiff
D(f, x, p, dx) = ForwardDiff.derivative(t->f(x .+ t .* dx, p), 0.)
d1F(x,p,dx1) = D((z, p0) -> F(z, p0), x, p, dx1)
d2F(x,p,dx1,dx2) = D((z, p0) -> d1F(z, p0, dx1), x, p, dx2)
d3F(x,p,dx1,dx2,dx3) = D((z, p0) -> d2F(z, p0, dx1, dx2), x, p, dx3)
jet = (F, J, d2F, d3F)
hopfpt = computeNormalForm(jet..., br, 3)
# newton parameters
optn_po = NewtonPar(verbose = true, tol = 1e-8, maxIter = 25)
# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.3, ds= -0.0001, dsmin = 1e-4, pMax = p_span[2], pMin=p_span[1], maxSteps = 1000, newtonOptions = (@set optn_po.tol = 1e-8), nev = 2, precisionStability = 4e-2, detectBifurcation = 0, plotEveryStep = 20, saveSolEveryStep=1)
Mt = 200 # number of sections
br_po, utrap = continuation(
jet..., br, 3, opts_po_cont,
PeriodicOrbitTrapProblem(M = Mt);
# specify the parameters for the initial guess. Otherwise they are extracted from opts_po_cont
ampfactor = 1., δp = 0.001,
# update the section every 1 steps. Very important
updateSectionEveryStep = 10,
# tangent predictor
tangentAlgo = BorderedPred(),
# how the jacobian of the Periodic orbit problem is computed. Here an optimised version for ODEs
linearPO = :Dense,
verbosity = 3, plot = true,
printSolution = (x, p) -> (xtt=reshape(x[1:end-1],2,Mt); return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = x[end])),
# function to plot the solution along the branch
plotSolution = (x, p; k...) -> begin
xtt = BK.getTrajectory(p.prob, x, p.p)
plot!(xtt.t, xtt.u[1,:]; label = "u1", markersize = 2, k...)
plot!(xtt.t, xtt.u[2,:]; label = "u2", k...)
plot!(br, br,subplot=1, putbifptlegend = false)
end,
# norm used for the computation
normC = x->norm(x,Inf)
)
plot(br, br_po, markersize = 3)
plot!(br_po.param, br_po.min, label = "")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment