비선형 ODE 시스템의 분기 플롯을 만드는 방법

Aug 21 2020

나는 박테리아 성장을 모델링하고 있습니다. 나는 ODE의 매우 어려운 비선형 시스템을 가지고 있습니다. 나는 경쟁으로 모델링하고있는 2 개의 박테리아를 가지고 있습니다. 따라서 모든 박테리아에 대해 2 개의 ODE가 있습니다. 왜냐하면 움직일 때 다른 것이되기 때문입니다 (P en S 방정식). 이제 매개 변수 분석을 시도하고 있으며 x 축에 매개 변수 중 하나 (a0, b0, d1, k, awmin, tmin)를 사용하여 분기 플롯을 만들려고 며칠 동안 노력해 왔습니다. 그러나 지금까지 시도한 것은 아무것도 작동하지 않았습니다. 좀 도와 주시고 어떻게 작동하는지 보여 주 시겠어요? 여기에 박테리아의 코드가 있습니다. 분기는 t = 150 부근에 자리를 잡아야합니다. 저는 이것에 대해 많은 문제를 겪었습니다.

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

답변

2 AlexTrounev Aug 23 2020 at 14:47

Module다음과 같이 파라 메트릭 연구에 사용할 수 있습니다.

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]

플롯하기 위해 Solution우리는 다음 f과 같이 사용 합니다.

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

Solution주어진 시간 동안 매개 변수의 함수로 플롯하려면 t=150새 함수를 정의합니다.

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