4
$\begingroup$

I have a differential equation system as follow:

B01 = {P -> 70, R -> 5, W -> 200, E1 -> 500, C1 -> 300, S -> 1, H -> 800};

G11 = (E1 - R*P - C1)*S /. B01
G12 = (E1 - H)*S /. B01
G13 = -C1*S /. B01
G14 = -H*S /. B01

C11 = (R*P - W)*S /. B01
C12 = -W*S /. B01
C13 = 0;
C14 = 0;

UG1[t_] := y[t] (G11) + (1 - y[t]) (G13)
UG2[t_] := y[t] (G12) + (1 - y[t]) (G14)

UC1[t_] := x[t] (C11) + (1 - x[t]) (C12)
UC2[t_] := x[t] (C13) + (1 - x[t]) (C14)

Solution1 = NDSolve[{x'[t] == x[t] (1 - x[t]) (UG1[t] - UG2[t]), y'[t] ==  y[t] (1 - y[t]) (UC1[t] - UC2[t]), y[0] == 0.25, 
    x[0] == 0.35}, {x, y}, {t, 0, 0.1}] /. B01

How I can plot:

  1. x[t] for x[0] in range of [0, 1] with interval 0.1 (for example: x[0] = {0, 0.1, 0.2, 03, ...})

  2. y[t] for y[0] in range [0, 1] with interval 0.1

  3. and at the end a parametric plot of $x$ vs $y$ similar to the picture below:

Pictures of 1 , 2 , 3

$\endgroup$

3 Answers 3

6
$\begingroup$
Clear["Global`*"]

B01 = {P -> 70, R -> 5, W -> 200, E1 -> 500, C1 -> 300, S -> 1, 
   H -> 800};
G11 = (E1 - R*P - C1)*S;
G12 = (E1 - H)*S;
G13 = -C1*S;
G14 = -H*S;
C11 = (R*P - W)*S;
C12 = -W*S;
C13 = 0;
C14 = 0;
UG1[t_] := y[t] (G11) + (1 - y[t]) (G13)
UG2[t_] := y[t] (G12) + (1 - y[t]) (G14)
UC1[t_] := x[t] (C11) + (1 - x[t]) (C12)
UC2[t_] := x[t] (C13) + (1 - x[t]) (C14)

eqns = {x'[t] == x[t] (1 - x[t]) (UG1[t] - UG2[t]), 
     y'[t] == y[t] (1 - y[t]) (UC1[t] - UC2[t]), y[0] == y0, 
     x[0] == x0} /. B01 // Simplify;

Solution1 = ParametricNDSolve[eqns, {x, y},
   {t, 0, 1/10}, {x0, y0}];

Plot[Evaluate[Flatten[Table[x[x0, y0][t],
     {x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}], 1] /. Solution1],
 {t, 0, 1/10},
 PlotRange -> All,
 Frame -> True,
 FrameLabel -> (Style[#, 14] & /@ {t, x}),
 AspectRatio -> 1]

enter image description here

Plot[Evaluate[Flatten[Table[y[x0, y0][t],
     {x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}], 1] /. Solution1],
 {t, 0, 1/10},
 PlotRange -> All,
 Frame -> True,
 FrameLabel -> (Style[#, 14] & /@ {t, y}),
 AspectRatio -> 1]

enter image description here

ParametricPlot[Evaluate[Flatten[Table[{x[x0, y0][t], y[x0, y0][t]},
     {x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}], 1] /. Solution1],
 {t, 0, 1/10},
 PlotRange -> All,
 Frame -> True,
 FrameLabel -> (Style[#, 14] & /@ {x, y}),
 AspectRatio -> 1,
 PlotPoints -> 75,
 MaxRecursion -> 5]

enter image description here

EDIT: Coordinating the PlotStyles

styles = ColorData[97] /@ Range[11];

Plot[
 Evaluate[
  Flatten[
    Table[x[x0, y0][t],
     {y0, 0, 1, 1/10}, {x0, 0, 1, 1/10}],
    1] /. Solution1],
 {t, 0, 1/10},
 PlotStyle -> styles,
 PlotRange -> All,
 Frame -> True,
 FrameLabel -> (Style[#, 14] & /@ {t, x}),
 AspectRatio -> 1]

enter image description here

Plot[
 Evaluate[
  Flatten[
    Table[y[x0, y0][t],
     {x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}],
    1] /. Solution1],
 {t, 0, 1/10},
 PlotStyle -> styles,
 PlotRange -> All,
 Frame -> True,
 FrameLabel -> (Style[#, 14] & /@ {t, y}),
 AspectRatio -> 1]

enter image description here

ParametricPlot[
 Evaluate[
  Flatten[
    Table[{x[x0, y0][t], y[x0, y0][t]},
     {x0, 0, 1, 1/10}, {y0, 0, 1, 1/10}],
    1] /. Solution1],
 {t, 0, 1/10},
 PlotStyle -> styles,
 PlotRange -> All,
 Frame -> True,
 FrameLabel -> (Style[#, 14] & /@ {x, y}),
 AspectRatio -> 1,
 PlotPoints -> 75,
 MaxRecursion -> 5]

enter image description here

$\endgroup$
1
  • $\begingroup$ Thanks Bob, This is what I need. $\endgroup$
    – Ahmad
    Jul 7 at 5:29
6
$\begingroup$

A quick way to get #3 is using StreamPlot:

StreamPlot[{x[t] (1 - x[t]) (UG1[t] - UG2[t]), 
  y[t] (1 - y[t]) (UC1[t] - UC2[t])}, {x[t], 0, 1}, {y[t], 0, 1}]

enter image description here

$\endgroup$
1
  • $\begingroup$ It is amazing, thanks Chris. $\endgroup$
    – Ahmad
    Jul 7 at 5:30
4
$\begingroup$

When I saw this question, I first thought StreamPlot, but it doesn't do the other graphs. Then I though wouldn't it be cool to convert an NDSolve solution to BezierCurve. I thought it would be cool because I had never done it. But after doing that, which might become a separate Q&A, I thought that the following would be cooler, even though it's not as robust ParametricNDSolveValue (see @BobHanlon's answer).

OP's setup in lowercase:

b01 = {p -> 70, r -> 5, w -> 200, e1 -> 500, c1 -> 300, s -> 1, 
   h -> 800};

g11 = (e1 - r*p - c1)*s /. b01;
g12 = (e1 - h)*s /. b01;
g13 = -c1*s /. b01;
g14 = -h*s /. b01;

c11 = (r*p - w)*s /. b01;
c12 = -w*s /. b01;
c13 = 0;
c14 = 0;

ug1[t_] := y[t] (g11) + (1 - y[t]) (g13);
ug2[t_] := y[t] (g12) + (1 - y[t]) (g14);

uc1[t_] := x[t] (c11) + (1 - x[t]) (c12);
uc2[t_] := x[t] (c13) + (1 - x[t]) (c14);

ics =(* we'll use the mesh coordinates as initial conditions *)
 DiscretizeRegion[Disk[{1/2, 1/2}, 1/2], 
  MaxCellMeasure -> 0.01]; solution2 = NDSolveValue[
  {x'[t] == x[t] (1 - x[t]) (ug1[t] - ug2[t]),
   y'[t] == y[t] (1 - y[t]) (uc1[t] - uc2[t]),
   (* we'll use the mesh coordinates all at once!: *)
   {x[0], y[0]} == Transpose@MeshCoordinates@ics},
  {x, y}, {t, -0.07, 0.07}];
(* dimensions of Through[solution2@"ValuesOnGrid"]
 * are {coord, step, ic} = 2 x 250 x 88 *)
systraj = Transpose[ (* transpose dims to {ic, step, coord} *)
    Through[solution2@"ValuesOnGrid"], {3, 2, 1}] //
   PadLeft[ #,  (* add t coord to x,y coords *)
     {Automatic, Automatic, 3}, (* {ic, step, new coords} *)
     solution2[[1]]@"Grid" (* use t from x sol, = same t as y *)
     ] &;
Dimensions@systraj
(*  {88, 250, 3}  *)

Graphics:

We can project the 3D trajectories onto 2D planes:

Table[
 Graphics[{
   Riffle[
    colors,
    Line /@ systraj[[All, All, parts]]]
   }, Options@Plot],
 {parts, {{1, 2}, {1, 3}, {2, 3}}}]

enter image description here

Or we can project the trajectories onto planes in 3D:

proj // ClearAll;
proj[k_, offsets_ : {-0.1, -0.3, -0.3}] :=(* 
  project onto coordinate plane *)
  TranslationTransform[
    ReplacePart[{0., 0., 0.}, k -> offsets[[k]]]] . 
   ScalingTransform[ReplacePart[{1., 1., 1.}, k -> 0.]];

colors = Hue[ (* VertexColors for the ICs *)
   Rescale[ArcTan @@ #, {-Pi, Pi}, {0, 1}],
   2 Norm@#,
   1
   ] &@(# - {1/2, 1/2}) & /@ MeshCoordinates@ics;

plall = Graphics3D[{
   Opacity[0.5], Thickness@0.003,
   Riffle[
    colors,
    Line /@ systraj],
   Opacity[0.2],
   Table[Riffle[
     colors,
     Line /@ proj[k]@systraj], {k, 3}],
   Append[ (* Initial conditions Disk[] *)
    ReplacePart[
        #,
        1 -> PadLeft[First@#, {Automatic, 3}, 0.]
        ] &@First@Show@ics /. {_Directive :> EdgeForm[Gray]},
    VertexColors -> colors]
   },
  BoxRatios -> {1, 1, 1},
  Axes -> True,
  AxesLabel -> {t, x, y},
  ViewPoint -> {2, 2.3, 1.8}, ViewVertical -> {0, 0, 1}]

Mathematica graphics

Like the plots in the OP:

pl3d = Graphics3D[{
    Opacity[0.5], Thickness@0.003,
    Riffle[
     colors,
     Line /@ systraj],
    Opacity[0.2],
    Append[(* Initial conditions Disk[] *)
     ReplacePart[
         #,
         1 -> PadLeft[First@#, {Automatic, 3}, 0.]
         ] &@First@Show@ics /. {_Directive :> EdgeForm[Gray]},
     VertexColors -> colors]
    },
   BoxRatios -> {1, 1, 1},
   Axes -> True,
   AxesLabel -> {t, x, y},
   ViewPoint -> {2, 2.3, 1.8}, ViewVertical -> {0, 0, 1}];

Grid[
   {{Graphics[Inset@#2, ImageSize -> 170], 
     Graphics[Inset@#1, ImageSize -> 340], 
     SpanFromLeft}, {Graphics[Inset@#3, ImageSize -> 170], 
     SpanFromAbove, SpanFromBoth}},
   Frame -> All, Alignment -> Top, Spacings -> {0, 0}
   ] & @@
 Table[Show[pl3d, opts],
  {opts, Transpose@{
     Thread[ViewPoint -> DiagonalMatrix[(-Infinity)^Range[2, 4]]], 
     Thread[PlotRange -> 
         {All,  {{0, 0.07}, All, All},  {{0, 0.07}, All, All}}],
     Thread[AxesEdge -> {
        { None,   {1, -1}, {1, -1}},
        {{1, -1},  None,   {-1, 1}},
        {{-1, 1}, {-1, 1} , None  }}],
     Thread[ImageSize -> {400, 200, 200}]}}]

Mathematica graphics

Without the projections, the graphics are smaller and easier to manipulate. If one zooms in on the initial condition disk, one can see the dynamics of the system better:

Show[pl3d, PlotRange -> {0.02 {-1, 1}, All, All}]
$\endgroup$
6
  • $\begingroup$ Whoa! I especially like that last one. $\endgroup$
    – Chris K
    Jul 7 at 23:52
  • $\begingroup$ oh my God, I am very impressed with your solutions. Your grapghs are very interesting and amazing. I appreciate and thank you very much for sharing these answers. $\endgroup$
    – Ahmad
    Jul 8 at 15:00
  • $\begingroup$ I thank you again for your solutions and plots Michael, Actually in my model time is time and it cant be negative, so we cant have negative parts of time in plots. $\endgroup$
    – Ahmad
    Jul 8 at 15:17
  • $\begingroup$ @Ahmad Thanks & you're welcome. As for time: (1) Since your model is autonomous, time is arbitrary. Solutions are invariant under translation $t\mapsto t+t_0$. So any solution may be translated to a time interval that is nonnegative. (2) You can always change the interval back to {t, 0, 0.1}. Note in the second-to-last graphic, the x and y vs. t plots are for nonnegative t. Using a symmetric interval was the easiest way to make the x vs. y plot to show Chris K's StreamPlot; using negative t does not invalidate the plot because the model is autonomous. (3)... $\endgroup$
    – Michael E2
    Jul 8 at 15:58
  • $\begingroup$ (3) My principal goal was to show some techniques that you or another user might adapt. And to make some pretty pictures. :) $\endgroup$
    – Michael E2
    Jul 8 at 15:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.