Come creare un diagramma di biforcazione di un sistema non lineare di ODE

Question

Sto modellando la crescita batterica. Ho un sistema non lineare abbastanza difficile di ODE. Ho 2 batteri che sto modellando con la concorrenza, oltre a trasferirli in un posto diverso. quindi ci sono 2 ODE per ogni batterio perché quando si muove sarà uno diverso (equazione P en S). Ora sto cercando di fare un'analisi dei parametri e ho provato per giorni a creare un grafico di biforcazione con sull'asse x uno dei parametri (a0, b0, d1, k, awmin, tmin). Ma niente di quello che ho provato finora ha funzionato. Ragazzi, potete aiutarmi un po 'e mostrarmi come funziona? ecco il codice dei batteri, la biforcazione dovrebbe avere un punto intorno a t = 150, ho avuto molti problemi con questo.

sp1 = D[G[t], t] == 
   Growth*G[t]*((Sp - G[t] - a0*B[t] )/Sp) + P - S  ;  
sp2 = D[M[t], 
    t] == (Growth2 - Kill)*M[t]*((Ss - M[t] - a0*A[t] )/Ss) + S - 
    P   ; 
P := M[t]* d1 ;  
S := G[t]*s*((G[t] + a0*B[t])/Sp) + G[t]*k     ;  
sp3 = D[B[t], t] == Growth*B[t]*((Sp - B[t] - b0*G[t] )/Sp) + P2 - S2;
sp4 = D[A[t], t] == 
   Growth2* A[t]*((Ss - A[t] - b0*M[t] )/Ss) + S2 - P2;
P2 := A[t]* d1;       
S2 := B[t]*0.01*((B[t] + b0*G[t] )/Sp) + B[t]*k;    
(*Parameters*)

Growth = 14.8*((T - tmin)*(Aw - Amin))^2  ;  
Growth2 = 7.4*((T - tmin)*(Aw - Amin))^2  ;
T := 4 + 18.09*Sin[0.01016*t + 0.3418];
Aw := 0.97 + 0.0001*t   ;
Kill = 1.5;
tmin = 12;
Amin = 0.95;
Sp = 1000;
Ss = 500;
a0 = 0.4;
b0 = 0.21;
d1 = 0.05;
k = 0.01;
s = 0.001;
all = { sp1, sp2, sp3, sp4};  
init1 = {G[0] == 60, M[0] == 3, B[0] == 50, A[0] == 30}; 
Solution = 
  NDSolveValue[{all, init1}, {G[t], M[t], B[t], A[t]}, {t, 0, 200}];  
Plot[Solution, {t, 0, 200}, PlotStyle -> {Red, Pink, Blue, Cyan}] 
ParametricPlot[{Solution[[1]], Solution[[3]]}, {t, 0, 200}, 
 PlotRange -> {{0, 100}, {0, 100}}, PlotStyle -> Red, 
 AspectRatio -> 1, PlotLabel -> "Phase plot", AxesLabel -> {"G", "B"}]

Answer 1

Possiamo usare Moduleper la ricerca parametrica come segue

f[a0p_, b0p_, d1p_, kp_, Ap_, tp_] := 
 Module[{a0 = a0p, b0 = b0p, d1 = d1p, k = kp, Amin = Ap, tmin = tp}, 
  Kill = 1.5;
  (*tmin=12;
  Amin=0.95;*)
  Sp = 1000;
  Ss = 500;
  (*a0=0.4;
  b0=0.21;
  d1=0.05;
  k=0.01;*)
  s = 0.001; 
  sp1 = D[G[t], t] == Growth*G[t]*((Sp - G[t] - a0*B[t])/Sp) + P - S;
  sp2 = D[M[t], 
     t] == (Growth2 - Kill)*M[t]*((Ss - M[t] - a0*A[t])/Ss) + S - P;
  P := M[t]*d1;
  S := G[t]*s*((G[t] + a0*B[t])/Sp) + G[t]*k;
  sp3 = D[B[t], t] == Growth*B[t]*((Sp - B[t] - b0*G[t])/Sp) + P2 - S2;
  sp4 = D[A[t], t] == 
    Growth2*A[t]*((Ss - A[t] - b0*M[t])/Ss) + S2 - P2;
  P2 := A[t]*d1;
  S2 := B[t]*0.01*((B[t] + b0*G[t])/Sp) + B[t]*k; 
  Growth = 14.8*((T - tmin)*(Aw - Amin))^2;
  Growth2 = 7.4*((T - tmin)*(Aw - Amin))^2;
  T := 4 + 18.09*Sin[0.01016*t + 0.3418];
  Aw := 0.97 + 0.0001*t; all = {sp1, sp2, sp3, sp4};
  init1 = {G[0] == 60, M[0] == 3, B[0] == 50, A[0] == 30};
  Solution = 
   NDSolveValue[{all, init1}, {G[t], M[t], B[t], A[t]}, {t, 0, 200}]; 
  Solution]

Per tracciare Solutionusiamo fcome

sol = f[.4, .21, .05, .01, .95, 12];
{Plot[sol, {t, 0, 200}, PlotStyle -> {Red, Pink, Blue, Cyan}],
 ParametricPlot[{sol[[1]], sol[[3]]}, {t, 0, 200}, 
  PlotRange -> {{0, 100}, {0, 100}}, PlotStyle -> Red, 
  AspectRatio -> 1, PlotLabel -> "Phase plot", 
  AxesLabel -> {"G", "B"}]}

Per tracciare Solutionin funzione dei parametri per un dato tempo t=150definiamo una nuova funzione, ad esempio,

solp[ap_] := f[ap, .21, .05, .01, .95, 12] /. t -> 150;
solp1[bp_] := f[.4, bp, .05, .01, .95, 12] /. t -> 150;
var = {G, M, B, A}; Table[
 Plot[solp[ap][[i]], {ap, .1, .5}, Frame -> True, 
  FrameLabel -> {"a0", var[[i]]}], {i, 4}]


Table[Plot[solp1[bp][[i]], {bp, .1, .5}, Frame -> True, 
  FrameLabel -> {"b0", var[[i]]}], {i, 4}]