How to make program for phase portrait?

Question

I need to add a phase portrait to my program. The program is below, and I need to create a graph like the one in the picture. My step size needs to be 5, from x0,y0 = 10 to x0,y0 = 50.

a := 0.1;
b := 0.2;
c := 0.1;

K := 80;
L := 70;

f[x_, y_] := a*x[t]*(1 - x[t]/K);
g[x_, y_] := b*y[t]*(1 - y[t]/(L - c*x[t]));

x0 := 10;
y0 := 10;
tmax := 100;


rjesenje :=
 NDSolve[{x'[t] == f[x, y], y'[t] == g[x, y], x[0] == x0,
   y[0] == y0}, {x, y}, {t, 0, tmax}]

ParametricPlot[Evaluate[{x[t], y[t]} /. rjesenje], {t, 0, tmax},
 PlotRange -> All, AxesOrigin -> {0, 0},
 AxesLabel -> {gustoća populacije x, gustoća populacije y}]

velicinax := x[t] /. Flatten[rjesenje][[1]]
velicinay := y[t] /. Flatten[rjesenje][[2]]

Plot[Evaluate[{velicinax, velicinay}], {t, 0, tmax}, PlotRange -> All,
  AxesOrigin -> {0, 0}, AxesLabel -> {vrijeme, gustoća populacije xy}]

How can I create a phase portrait like this?

Answer 1

You may do as follows. The function below shows a single trajectory:

traj1[eq1_, eq2_, point_, col_, tmax_, n_] :=
        Module[{p0, p1, p2, p3, tau, s},
          eq3 = x[0] == point[[1]]; 
             eq4 = y[0] == point[[2]];
                s = NDSolve[{eq1, eq2, eq3, eq4}, {x, y}, {t, tmax}, 
           Method -> "StiffnessSwitching"];
     tau = tmax/n;
p0 = ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, tmax}, 
 PlotStyle -> col];
    p1 = Graphics[{col,
        Arrow[{Evaluate[{x[tau], y[tau]} /. s] // Flatten, 
             Evaluate[{x[tau + 0.1], y[tau + 0.1]} /. s] // Flatten}]}];
          p2 = Graphics[{col,
    Arrow[{Evaluate[{x[2 tau], y[2 tau]} /. s] // Flatten, 
Evaluate[{x[2 tau + 0.1], y[2 tau + 0.1]} /. s] // Flatten}]}];
            p3 = Graphics[{col,                  
  Arrow[{Evaluate[{x[3 tau], y[3 tau]} /. s] // Flatten, 
    Evaluate[{x[3 tau + 0.1], y[3 tau + 0.1]} /. s] // Flatten}]}];
             Show[{p0, p1, p2, p3}] ];

Arguments eq1, and eq2 are equations, point is a list {x0,y0} giving the initial point for a given trajectory, col is the trajectory color, tmax is the same as in your code and n indirectly controls the number of the arrowheads per trajectory, the value between 5 and 10 is recommended.

Now, this is your input:

a = 0.1;
b = 0.2;
c = 0.1;
k = 80;
L = 70;

eq1 = x'[t] == a*x[t]*(1 - x[t]/k);
eq2 = y'[t] == b*y[t]*(1 - y[t]/(L - c*x[t]));

Take care that I use k, rather than K, which is important. Further this builds the phase portrait with four trajectories:

  Manipulate[
 Show[{traj1[eq1, eq2, {20, 10}, Red, tmax, 10], 
   traj1[eq1, eq2, {0, 10}, Blue, tmax, 10], 
   traj1[eq1, eq2, {20, 20}, Green, tmax, 10], 
   traj1[eq1, eq2, {5, 20}, Purple, tmax, 10]}, PlotRange -> All],
 {tmax, 10, 100}]

You might want to add more trajectories just by adding more traj1 functions under the Show statement. Manipulation is added to see the effect of different tmax. You should obtain this on the screen:

This is as the first, simple shot. One can, of course, think also about something more complex and flexible, if needed. Have fun!

Later edit: To address your question, this is a list of the initial conditions:

 lst = Flatten[Table[{10 + 20 i, 10 + 20 j}, {i, 0, 5}, {j, 0, 5}], 1];

Note that I did not follow your request:" of values 10 to 50 with step 5" for several reasons: a) the trajectories are much too dense in this case and one sees nothing, b) the right part of the phase portrait is not visible.
Then (provided you have already evaluated all previous statements) this:

 Show[traj1[eq1, eq2, #, Hue[N[Norm[#]/110]], 100, 7] & /@ lst, 
 PlotRange -> {{0, 100}, {0, 100}}]

gives you the following phase portrait:

The trajectory color is controlled by the function Hue and depends upon the absolute value of the radius-vector of the initial point: Sqrt(x0^2+y0^2). You might want to play with Hue to achieve the color scheme you desire (and which I did not really understand).

There may be also an approach based on a random set of initial points as follows:

lst2 = RandomReal[{0, 150}, {35, 2}];

Show[traj1[eq1, eq2, #, Hue[N[#[[1]]/Max[Transpose[lst2][[1]]]]], 100,
 7] & /@ lst2, PlotRange -> {{0, 100}, {0, 100}}]

which returns this:

where the colors depend upon x0, but independent of y0. There are plenty of other possibilities, including approaches of others, as pointed out by Michael E2 in the discussion.

So, generally, have fun!