4

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 - RP - C1)S /. B01 G12 = (E1 - H)S /. B01 G13 = -C1S /. B01 G14 = -H*S /. B01

C11 = (RP - 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

MarcoB
  • 67,153
  • 18
  • 91
  • 189
Ahmad
  • 77
  • 5

3 Answers3

6
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

Bob Hanlon
  • 157,611
  • 7
  • 77
  • 198
6

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

Chris K
  • 20,207
  • 3
  • 39
  • 74
4

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 - rp - c1)s /. b01; g12 = (e1 - h)s /. b01; g13 = -c1s /. b01; g14 = -h*s /. b01;

c11 = (rp - 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}]
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Whoa! I especially like that last one. – Chris K Jul 07 '22 at 23:52
  • 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. – Ahmad Jul 08 '22 at 15:00
  • 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. – Ahmad Jul 08 '22 at 15:17
  • @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)... – Michael E2 Jul 08 '22 at 15:58
  • (3) My principal goal was to show some techniques that you or another user might adapt. And to make some pretty pictures. :) – Michael E2 Jul 08 '22 at 15:58
  • I personally enjoyed your techniques and plots, all of them are very interesting for me, Thanks @Michael E2. – Ahmad Jul 08 '22 at 16:16