Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# rveltz/catalystexp1.jl Secret

Created May 12, 2021
 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 = "")
to join this conversation on GitHub. Already have an account? Sign in to comment