Cómo hacer una gráfica de bifurcación de un sistema no lineal de EDO

Question

Estoy modelando el crecimiento bacteriano. Tengo un sistema no lineal de EDO bastante difícil. Tengo 2 bacterias que estoy modelando con la competencia, además de transferirlas a un lugar diferente. por lo que hay 2 EDO para cada bacteria porque cuando se mueva será una diferente (ecuación P en S). Ahora estoy tratando de hacer un análisis de parámetros y he estado intentando durante días hacer una gráfica de bifurcación con uno de los parámetros en el eje x (a0, b0, d1, k, awmin, tmin). Pero nada de lo que he intentado hasta ahora funcionó. ¿Pueden ayudarme un poco y mostrarme cómo funciona? aquí está el código de la bacteria, la bifurcación debería tener un lugar alrededor de t = 150, he tenido muchos problemas con esto.

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

Podemos utilizar Modulepara la investigación paramétrica de la siguiente manera

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]

Para trazar Solutionusamos fcomo

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"}]}

Para graficar Solutioncomo una función de parámetros para un tiempo dado t=150, definimos una nueva función, por ejemplo,

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}]