Mathematica做2物种竞争以及3物种竞争的动力系统仿真

目录

• upd 2021-02-06 完成第一版改造
  • 3物种
  • 2物种 魔改前
  • 2物种 魔改后
• 参考

upd 2021-02-06 完成第一版改造

魔改Wolfram Demonstration的代码
建模利用的是这个Competitive Lotka–Volterra equations
本篇博客园文章没有对物种和种群进行区分

\[\frac{d x_{i}}{d t}=r_{i} x_{i}\left(1-\frac{\sum_{j=1}^{N} \alpha_{i j} x_{j}}{K_{i}}\right) \]

自治系统理解为微分方程不显含时间t

3物种

 Manipulate[

sol=NDSolve[

{

x1'[t]==r1*x1[t]*(1-\[Alpha]11*x1[t]-\[Alpha]12 x2[t]-\[Alpha]13 x3[t]),

x2'[t]==r2*x2[t]*(1-\[Alpha]21*x1[t]-\[Alpha]22 x2[t]-\[Alpha]23 x3[t]),

x3'[t]==r3*x3[t]*(1-\[Alpha]31*x1[t]-\[Alpha]32 x2[t]-\[Alpha]33 x3[t]),

x1[0]==x10,

x2[0]==x20,

x3[0]==x30

},

{x1[t],x2[t],x3[t]},

{t,0,time}

];

Pane[

Grid[{{

Show[

{

ParametricPlot3D[{x1[t],x2[t],x3[t]}/.sol,{t,0,time},

PlotStyle->{Thick,Blue}],

Graphics3D[

{Red,PointSize[Large],Point[{x1[t],x2[t],x3[t]}/.sol/.{t->0}],Green,Point[{x1[t],x2[t],x3[t]}/.sol/.{t->time}]}

]

},

If[range==="fixed",PlotRange->{{0,x1max},{0,x2max},{0,x3max}},PlotRange->All],

BoxRatios->{1,1,1},

If[values,Ticks->Automatic,Ticks->None],

AxesLabel->(Style[#,Blue]&/@(Row/@Transpose[{Table[Subscript[Style["X",Italic],i],{i,3}],If[label,{"\n(Lowest\ntrophic level)","\n(Intermediate\ntrophic level)","\n(Top\npredators)"},{"","",""}]}])),

ImageSize->{400,400},ImagePadding->35

]},

{Text@Grid[{

{

Grid[{

{Row[{

Subscript[Style["X",Italic], 3],

" ("<>ToString[Round[100x3[t]/(x1[t]+x2[t]+x3[t])/.sol/.{t->time}][[1]]]<>"%)"

}]},

{Row[{

Subscript[Style["X",Italic], 2],

" ("<>ToString[Round[100x2[t]/(x1[t]+x2[t]+x3[t])/.sol/.{t->time}][[1]]]<>"%)"

}]},

{Row[{

Subscript[Style["X",Italic], 1],

" ("<>ToString[Round[100x1[t]/(x1[t]+x2[t]+x3[t])/.sol/.{t->time}][[1]]]<>"%)"

}]}

}],



BarChart[Round@Flatten[{x1[t],x2[t],x3[t]}/.sol/.{t->time}],

BarSpacing->0,

BarOrigin->Left,

ChartStyle->{{RGBColor[.1,.9,.1],RGBColor[.4,.6,.4],RGBColor[.7,.3,.7]}},

Axes->None,

LabelingFunction->Left,

AspectRatio->.25,

ImageSize->{200,100}

]

}



},

Alignment->{Right,Left}]}



}],

ImageSize->400

],

Style["initial variable values",Bold],

{{x10,2000,Subscript["X",1]},1,7000,.01,Appearance->"Labeled",ImageSize->Tiny},

{{x20,175,Subscript["X",2]},1,7000,.01,Appearance->"Labeled",ImageSize->Tiny},

{{x30,200,Subscript["X",3]},1,7000,.01,Appearance->"Labeled",ImageSize->Tiny},



Style["parameter values",Bold],

{{r1,2.84,Subscript["r",1]},0,5,.01,Appearance->"Labeled",ImageSize->Tiny},

{{r2,1.5,Subscript["r",2]},0,5,.01,Appearance->"Labeled",ImageSize->Tiny},

{{r3,.62,Subscript["r",3]},0,5,.01,Appearance->"Labeled",ImageSize->Tiny},

{{\[Alpha]11,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{\[Alpha]12,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{\[Alpha]13,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{\[Alpha]21,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{\[Alpha]22,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{\[Alpha]23,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{\[Alpha]31,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{\[Alpha]32,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},
{{\[Alpha]33,.22},0,1,.01,Appearance->"Labeled",ImageSize->Tiny},



Style["time period",Bold],

{{time,10,""},1,100,.01,Appearance->"Labeled",ImageSize->Tiny},



Style["maximum values on axes",Bold],

{{range,"floating",""},{"fixed","floating"}},

{{x1max,3100,Subscript["X",1]},0,20000,.01,Appearance->"Labeled",ImageSize->Tiny},

{{x2max,300,Subscript["X",2]},0,5000,.01,Appearance->"Labeled",ImageSize->Tiny},

{{x3max,500,Subscript["X",3]},0,5000,.01,Appearance->"Labeled",ImageSize->Tiny},



Style["label trophic levels",Bold],

{{label,False,""},{True,False}},



Style["display values on axes",Bold],

{{values,True,""},{True,False}},

TrackedSymbols->True,SynchronousUpdating->True,

ControlPlacement->Left,AutorunSequencing->{2,5,9,13,15}

]

2物种 魔改前

Manipulate[
 Module[{pic, sol},
  pic = VectorPlot[{eps1*x - s1*x^2 - \[Alpha]1*y*x, 
     eps2*y - s2*y^2 - \[Alpha]2*y*x}, {x, 0, 4}, {y, 0, 4}, 
    StreamPoints -> {{p[[1]], p[[2]]}}, StreamStyle -> {Red, Thick}, 
    VectorScale -> {Tiny, Tiny, None}, VectorStyle -> "Segment", 
    ImageSize -> {350, 350}, 
    PlotLabel -> 
     Column[{Row[{Style["x", Italic], "'", " = ", eps1*x, " - ", 
         s1*x^2, " - ", \[Alpha]1*y*x}], 
       Row[{Style["y", Italic], "'", " = ", eps2*y, " - ", s2*y^2, 
         " - ", \[Alpha]2*y*x}], Spacer[12]}], 
    FrameLabel -> {Row[{"foxes ", Style["x", Italic]}], 
      Row[{"hawks ", Style["y", Italic]}]}];
  sol = NDSolve[{{x1'[t] == 
       eps1*x1[t] - s1*(x1[t])^2 - \[Alpha]1*y1[t]*x1[t], 
      y1'[t] == 
       eps2*y1[t] - s2*(y1[t])^2 - \[Alpha]2*y1[t]*x1[t]}, {x1[0] == 
       p[[1]], y1[0] == p[[2]]}}, {x1[t], y1[t]}, {t, 0, tend}];
  x2[t_] := x1[t] /. sol[[1, 1]];
  y2[t_] := y1[t] /. sol[[1, 2]];
  plotx = 
   Plot[x2[t], {t, 0, tend}, 
    PlotLabel -> Row[{"foxes ", Style["x", Italic]}], 
    PlotStyle -> Red, PlotRange -> {0, 10}, ImageSize -> {200, 100}, 
    AspectRatio -> .5];
  ploty = 
   Plot[y2[t], {t, 0, tend}, 
    PlotLabel -> Row[{"hawks ", Style["y", Italic]}], 
    PlotStyle -> Blue, PlotRange -> {0, 10}, ImageSize -> {200, 100}, 
    AspectRatio -> .5];
  Show[pic, ImageSize -> {370, 370}]],
 Style["parameters", Bold],
 {{eps1, 1, Subscript[\[Epsilon], 1]}, 0, 1, .1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 {{s1, .5, Subscript[\[Sigma], 1]}, 0, 1, .1, Appearance -> "Labeled",
   ImageSize -> Tiny},
 {{\[Alpha]1, .4, Subscript[\[Alpha], 1]}, 0, 1, .1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 "",
 {{eps2, 1, Subscript[\[Epsilon], 2]}, 0, 1, .1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 {{s2, .5, Subscript[\[Sigma], 2]}, 0, 1, .1, Appearance -> "Labeled",
   ImageSize -> Tiny},
 {{\[Alpha]2, .25, Subscript[\[Alpha], 2]}, 0, 1, .1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 {{plotx, Graphics[]}, ControlType -> None},
 {{ploty, Graphics[]}, ControlType -> None},
 "",
 Delimiter,
 {{tend, 10, "end time \!\(\*SubscriptBox[\(t\), \(max\)]\)"}, 10, 
  100, 5, Appearance -> "Labeled", ImageSize -> Tiny},
 Delimiter,
 Dynamic[Show[plotx], SynchronousUpdating -> False],
 Delimiter,
 Dynamic[Show[ploty], SynchronousUpdating -> False],
 {{p, {.5, .5}}, {0, 0}, {4, 4}, Locator},
 {{plotx, Graphics[]}, ControlType -> None},
 {{ploty, Graphics[]}, ControlType -> None},
 ControlPlacement -> Left,
 AutorunSequencing -> {1, 3},
 TrackedSymbols :> {eps1, eps2, s1, s2, \[Alpha]1, \[Alpha]2, tend, p},
 SynchronousUpdating -> False]

2物种 魔改后

Manipulate[
 Module[{pic, sol}, 
  pic = VectorPlot[{r1*x*(1 - \[Alpha]11*x - \[Alpha]12*y), 
     r2*y*(1 - \[Alpha]22*y - \[Alpha]21*x)}, {x, 0, 20}, {y, 0, 20}, 
    StreamPoints -> {{p[[1]], p[[2]]}}, StreamStyle -> {Red, Thick}, 
    VectorScale -> {Tiny, Tiny, None}, VectorStyle -> "Segment", 
    ImageSize -> {350, 350}, 
    PlotLabel -> 
     Column[{Row[{Style["x", Italic], "'", " = ", r1*x, "(", 1, 
         "-", \[Alpha]11*x, "-", \[Alpha]12*y, ")"}], 
       Row[{Style["y", Italic], "'", " = ", r2*y, "(", 1, 
         "-", \[Alpha]22*y, "-", \[Alpha]21*x, ")"}], Spacer[12]}], 
    FrameLabel -> {Row[{"Species1 ", Style["x", Italic]}], 
      Row[{"Species2 ", Style["y", Italic]}]}];
  sol = NDSolve[{{x1'[t] == 
       r1*x1[t]*(1 - \[Alpha]11*x1[t] - \[Alpha]12*y1[t]), 
      y1'[t] == 
       r2*y1[t]*(1 - \[Alpha]22*y1[t] - \[Alpha]21*x1[t])}, {x1[0] == 
       p[[1]], y1[0] == p[[2]]}}, {x1[t], y1[t]}, {t, 0, tend}];
  x2[t_] := x1[t] /. sol[[1, 1]];
  y2[t_] := y1[t] /. sol[[1, 2]];
  plotx = 
   Plot[x2[t], {t, 0, tend}, 
    PlotLabel -> Row[{"Species1 ", Style["x", Italic]}], 
    PlotStyle -> Red, PlotRange -> {0, 20}, ImageSize -> {200, 100}, 
    AspectRatio -> .5];
  ploty = 
   Plot[y2[t], {t, 0, tend}, 
    PlotLabel -> Row[{"Species2 ", Style["y", Italic]}], 
    PlotStyle -> Blue, PlotRange -> {0, 20}, ImageSize -> {200, 100}, 
    AspectRatio -> .5];
  Show[pic, ImageSize -> {370, 370}]], 
 Style["parameters", Bold], {{r1, 0.5, Subscript[r, 1]}, 0, 1, .1, 
  Appearance -> "Labeled", 
  ImageSize -> Tiny}, {{\[Alpha]11, 0.1, Subscript[\[Alpha], 11]}, 0, 
  1, .1, Appearance -> "Labeled", 
  ImageSize -> Tiny}, {{\[Alpha]12, 0.4, Subscript[\[Alpha], 12]}, 0, 
  1, .1, Appearance -> "Labeled", 
  ImageSize -> Tiny}, "", {{r2, 0.2, Subscript[r, 2]}, 0, 1, .1, 
  Appearance -> "Labeled", 
  ImageSize -> Tiny}, {{\[Alpha]22, 0.4, Subscript[\[Alpha], 22]}, 0, 
  1, .1, Appearance -> "Labeled", 
  ImageSize -> Tiny}, {{\[Alpha]21, 0.3, Subscript[\[Alpha], 21]}, 0, 
  1, .1, Appearance -> "Labeled", 
  ImageSize -> Tiny}, {{plotx, Graphics[]}, 
  ControlType -> None}, {{ploty, Graphics[]}, 
  ControlType -> None}, "", Delimiter, {{tend, 10, 
   "end time Subscript[t, max]"}, 10, 100, 5, Appearance -> "Labeled",
   ImageSize -> Tiny}, Delimiter, 
 Dynamic[Show[plotx], SynchronousUpdating -> False], Delimiter, 
 Dynamic[Show[ploty], 
  SynchronousUpdating -> False], {{p, {.5, .5}}, {0, 0}, {20, 20}, 
  Locator}, {{plotx, Graphics[]}, 
  ControlType -> None}, {{ploty, Graphics[]}, ControlType -> None}, 
 ControlPlacement -> Left, AutorunSequencing -> {1, 3}, 
 TrackedSymbols :> {r1, \[Alpha]11, \[Alpha]12, 
   r2, \[Alpha]22, \[Alpha]21, tend, p}, SynchronousUpdating -> False]

参考

https://demonstrations.wolfram.com/EcosystemDynamics/
https://demonstrations.wolfram.com/CompetingSpeciesModelWithTwoSpecies/