# How to make program for phase portrait?

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 == x0,
y == 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][]
velicinay := y[t] /. Flatten[rjesenje][]

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? • Please add your code in a copyable format, thereby making it easier for others to help you. You can use the {} button of the editor for proper formatting and check the editing help for additional informations. – Karsten 7. Jun 24 '15 at 11:04
• @Karsten7. is this okey? – ivona0909 Jun 24 '15 at 11:29
• @ivona0909 I have edited your question for readability, I hope this is what you meant to ask. – shrx Jun 24 '15 at 12:21
• Have you seen StreamPlot and EquationTrekker? Also, there are oodles of questions on phase portrait and phase diagram on site already. – Michael E2 Jun 24 '15 at 17:50
• @AlexeiBoulbitch This is great, but can it do it a a loop that would make the same as this in this picture; Trajectories are defined for different initial conditions x0 and y0 that range of values 10 to 50 with step 5. A specific color-coded the different initial values of y for the same the initial value of x. So, for example, dark blue marked the trajectory of x0 = 20, and y 0 = 10, 15, 20, 25, 30, 35, 40, 45 and 50th – ivona0909 Jun 24 '15 at 19:05

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 == point[];
eq4 = y == point[];
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.

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[#[]/Max[Transpose[lst2][]]]], 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!