💡💡💡 DYNAMO 2x2 💡💡💡
Phase Diagrams for Evolutionary Dynamics

User-defined Parameters

When you have finished specifying parameters, click on the rightmost vertical bar and press Enter (or Shift-Enter) to run the program.

Off[ General::spell1 ] ;  Off[ General::spell] ;

Choice of game

Choose a normal form game here.  To use one of the built-in payoff matrices, select the matrix on the right hand side of "U = ." (including the parentheses but not the semicolon); then the press the button corresponding to the game you desire.

({{12 Coordination , Pure Coordination, Matching Pennies , Prisoner' s Dilemma , Stag Hunt, Chain Store}})

U=(({{{1, -1}, {-1, 1}}, {{-1, 1}, {1, -1}}}));

Nonlinear population games

To choose a nonlinear game, set nonlineargame = 1 and define F1[x_,y_]  and F2[x_,y_] as you like. In doing so, the components of x should be entered as x[[1]], x[[2]],  and the components of y should be entered similarly.  The default value of F1[x_,y_] is A.y and F2[x_,y_] is B'.x, the payoff vector field for the normal form game U.

nonlineargame  = 0 ;

A=(({{U[[1, 1, 
                  1]], U[[1, 2, 
                  1]]}, {U[[2, 1, 
                  1]], U[[2, 2, 
                  1]]}})) ; B=(({{U[[1, 1, 
                  2]], U[[1, 2, 
                  2]]}, {U[[2, 1, 
                  2]], U[[2, 2, 
                  2]]}})) ;
F1[x_, y_]  := A . y ; 
F2[x_, y_]  := Transpose[B] . x ;  
    F[x_, y_] := {F1[x, y], F2[x, y]} ;

If the game you consider has a nontrivial component of Nash equilibria, the program will find a random subset of this component.  If you'd prefer not to compute the Nash equilibria, set findnashequilibria to 0.

findnashequilibria = 1 ;

Give the strategies one-character names to be printed on the phase diagram.

strategy[1] = "T" ; 
strategy[2] = "B" ; 
strategy[3] = "L" ;
strategy[4] = "R" ;

Payoff-related definitions

Clear[Fbar] ;  Clear[Fhat] ;   Clear[Fhatplus] ;

The population's average payoff

Fbar[1][x_, y_, 
          F1_] := (∑_(i = 1)^2) x[[
              i]] F1[x, y][[
              i]] ;

    
Fbar[2][x_, y_, 
          F2_] :=  (∑_(i = 1)^2) y[[
              i]] F2[x, y][[
              i]] ;

Excess payoff vector

Fhat[1][x_, y_, F1_] := 
        F1[x, y]  -  Fbar[1][x, y, F1] ;

    
Fhat[2][x_, y_, F2_] := 
        F2[x, y]  -  Fbar[2][x, y, F2] ;

Vector of positive parts of excess payoffs

Fhatplus[1][x_, y_, 
          F1_] :=  {Max[
            0, Fhat[1][x, y, 
                F1][[1]]], 
          Max[0, Fhat[1][x, y, 
                F1][[2]]]} ;

    
Fhatplus[2][x_, y_, 
          F2_] :=  {Max[
            0, Fhat[2][x, y, 
                F2][[1]]], 
          Max[0, Fhat[2][x, y, 
                F2][[2]]]} ;

Projected payoff vector

PhiF[x_, y_, F1_, F2_] := 
        Flatten[{F1[x, 
                y] - {1/2, 
                  1/2} (∑_(i = 1)^2) F1[x, y][[i]], 
            F2[x, y] - {1/2, 
                  1/2} (∑_(i = 1)^2) F2[x, y][[i]] }, 
          1];

Choice of dynamic

To specify a dynamic, select the rhs of "Dynamic = (.)" (excluding the semicolon), and then press the button corresponding to the dynamic you desire.  Definitions of the dynamics can be found in the closed group at the end of this section.
Details:  If you use Logit[η] or ILogit[η], be sure to enter a value for the noise parameter η in the closed group below.  Making  η very small generates close approximations of the best response dynamic.  If you use Projection and if you are not using Mathematica version 5 or above, do not take initial condtions on the boundary. EP (excess payoff dynamic), PC (pairwise comparison dynamic), ConvexComb, and Other must be specified in greater detail below.

Noise level for the logit and i-logit dynamics

η  =   .01 ;

({{Replicator, Logit[η], BNN, PD, Projection, Combined}, {MSReplicator, ILogit[η], EP, PC, , Other}})

Dynamic = Replicator ;

If the dynamic you consider has a nontrivial component of rest points, the program will compute and plot a random subset of this component.  If you set findrestpoints equal to 0, the computation and plotting is skipped.

findrestpoints = 1;

Definition of combination of two dynamics

To define a combination of two dynamics, specify the two dynamics to be combined, the weights  on the dynamics, and the characterization of the rest points of combined dynamic.  Nash is the default setting for the set or rest points.  This is correct, for example, with combinations of the replicator dynamic and an excess payoff dynamic like the BNN dynamic.  If the rest points aren't identical to the Nash equilibria for the dynamics you have specified, you can replace Nash with Automatic if the combined dynamic is smooth.  (If the combined dynamic isn't smooth, you're on your own!)  The last line of code formally defines the combined dynamic according to your specificcations.

dyn[1] = Replicator;
dyn[2] = BNN;
weight[1] =  .9;
weight[2] =  .1;
RPCharacterization[Combined] = Nash;
Combined[x_, y_, F1_, F2_] := 
        weight[1] dyn[1][x, y, F1, F2] + 
          weight[2] dyn[2][x, y, F1, F2];

Definitions of dynamics

Clear[Replicator] ;  
Clear[BNN]  ;
Clear[Logit] ;
Clear[rho] ;
Clear[PD] ; 
Clear[rho1] ;
Clear[SelMut] ;

Replicator Dynamic

Replicator[x_, y_, F1_, 
          F2_] := {x[[1]]  Fhat[
                  1][x, y, F1][[1]], 
          x[[2]]  Fhat[1][x, y, 
                F1][[2]], 
          y[[1]]  Fhat[2][x, y, 
                F2][[1]], 
          y[[2]]  Fhat[2][x, y, 
                F2][[2]]} ;

RPCharacterization[Replicator] = 
        Automatic ;

Maynard Smith Replicator Dynamic

MSReplicator[x_, y_, F1_, 
          F2_] := {(1/(Fbar[1]))[x, y, F1] 
            x[[1]]  Fhat[1][x, 
                y, F1][[1]], (1/(Fbar[1]))[x, y, F1] 
            x[[2]]  Fhat[1][x, 
                y, F1][[2]], (1/(Fbar[2]))[x, y, F2] 
            y[[1]]  Fhat[2][x, 
                y, F2][[1]], (1/(Fbar[2]))[x, y, F2] 
            y[[2]]  Fhat[2][x, 
                y, F2][[2]]} ;

RPCharacterization[MSReplicator] = 
        Automatic ;

Logit Dynamic

Logit[η_][x_, y_, F1_, 
          F2_] := {((E^(F1[x, y][[1]]/η))/((∑_(i = 1)^2) (E^(F1[x, y][[i]]/η)))) - 
            x[[1]], ((E^(F1[x, y][[2]]/η))/((∑_(i = 1)^2) (E^(F1[x, y][[i]]/η)))) - 
            x[[2]], ((E^(F2[x, y][[1]]/η))/((∑_(i = 1)^2) (E^(F2[x, y][[i]]/η)))) - 
            y[[1]], ((E^(F2[x, y][[2]]/η))/((∑_(i = 1)^2) (E^(F2[x, y][[i]]/η)))) - 
            y[[2]]} ;

RPCharacterization[Logit[_]] = 
        Automatic ;

I-Logit Dynamic

ILogit[η_][x_, y_, F1_, 
          F2_] := {((x[[1]] (E^(F1[x, y][[1]]/η)))/((∑_(i = 1)^2) x[[i]] (E^(F1[x, y][[i]]/η))))  - 
            x[[1]], ((x[[2]] (E^(F1[x, y][[2]]/η)))/((∑_(i = 1)^2) x[[i]] (E^(F1[x, y][[i]]/η)))) - 
            x[[2]], ((y[[1]] (E^(F2[x, y][[1]]/η)))/((∑_(i = 1)^2) y[[i]] (E^(F2[x, y][[i]]/η))))  - 
            y[[1]], ((y[[2]] (E^(F2[x, y][[2]]/η)))/((∑_(i = 1)^2) y[[i]] (E^(F2[x, y][[i]]/η)))) - 
            y[[2]]} ;

RPCharacterization[ILogit[_]] = 
        Automatic ;

BNN dynamic

BNN[x_, y_, F1_, 
          F2_] := {Fhatplus[1][x, y, 
                F1][[1]] - 
            x[[1]] (∑_(i = 1)^2) Fhatplus[1][x, y, F1][[
                  i]], Fhatplus[1][x, y, 
                F1][[2]] - 
            x[[2]] (∑_(i = 1)^2) Fhatplus[1][x, y, F1][[
                  i]], Fhatplus[2][x, y, 
                F2][[1]] - 
            y[[1]] (∑_(i = 1)^2) Fhatplus[2][x, y, F2][[
                  i]], Fhatplus[2][x, y, 
                F2][[2]] - 
            y[[2]] (∑_(i = 1)^2) Fhatplus[2][x, y, F2][[
                  i]]} ;

RPCharacterization[BNN] = Nash ;

Excess payoff dynamic

Here one specifies the vector field sigmatilde that defines an excess payoff dynamic.  This definition should be stated in terms of the excess payoff vector Fhat[i][x,y,Fi], which is defined as Fhat[i][x,y,Fi] = Fi[x,y] - Fbar[i][x,y,Fi], for each population i = 1, 2.  (Mathematica subtracts the scalar from the vector correctly.)  When computing rest points, the program assumes that sigmatilde is acute, so that the rest points of the dynamic are the Nash equilibria of the underlying game.

sigmatilde[x_, y_, F1_, 
          F2_] :=  {(Max[
                0, Fhat[1][x, y, 
                    F1][[1]]])^2, 
(Max[
                0, Fhat[1][x, y, 
                    F1][[2]]])^2, (Max[0, Fhat[2][x, y, 
                    F2][[1]]])^2, (Max[0, Fhat[2][x, y, 
                    F2][[2]]])^2};

EP[x_, y_, F1_, 
          F2_] := {sigmatilde[x, y, F1, 
                F2][[1]] - 
            x[[1]] (∑_(i = 1)^2) sigmatilde[x, y, F1, F2][[
                  i]], sigmatilde[x, y, F1, 
                F2][[2]] - 
            x[[2]] (∑_(i = 1)^2) sigmatilde[x, y, F1, F2][[
                  i]], sigmatilde[x, y, F1, 
                F2][[3]] - 
            y[[1]] (∑_(i = 3)^4) sigmatilde[x, y, F1, F2][[
                  i]], sigmatilde[x, y, F1, 
                F2][[4]] - 
            y[[2]] (∑_(i = 3)^4) sigmatilde[x, y, F1, F2][[
                  i]]} ;

RPCharacterization[EP] = Nash  ;

Pairwise Difference Dynamic

rho[1][x_, y_, 
          F1_] := {{0, 
            Max[F1[x, y][[2]] - F1[
                    x, y][[1]], 
              0]}, {Max[F1[x, 
                    y][[1]] - F1[x, 
                    y][[2]], 0], 
            0}} ;
rho[2][x_, y_, 
          F2_] := {{0, 
            Max[F2[x, y][[2]] - F2[
                    x, y][[1]], 
              0]}, {Max[F2[x, 
                    y][[1]] - F2[x, 
                    y][[2]], 0], 
            0}} ;

PD[x_, y_, F1_, F2_] := 
        Join[Flatten[
            Transpose[rho[1][x, y, F1]] .  
                x - {{x[[1]], 0}, {0, 
                    x[[2]]}} . rho[
                    1][x, y, F1] . {1, 1}], 
          Flatten[Transpose[rho[2][x, y, F2]] .  
                y - {{y[[1]], 0}, {0, 
                    y[[2]]}} . rho[
                    2][x, y, F2] . {1, 1}]] ;

RPCharacterization[PD] = Nash ;

Pairwise Comparison Dynamic

Here one specifies the vector field rho that defines an excess payoff dynamic.  When computing rest points, the program assumes that rho is separable and sign-preserving, so that the rest points of the dynamic are the Nash equilibria of the underlying game.

rho1[1][x_, y_, 
          F1_] := {{0, 
            Max[F1[x, y][[2]] - F1[
                    x, y][[1]], 
              0]}, {Max[F1[x, 
                    y][[1]] - F1[x, 
                    y][[2]], 0], 
            0}} ;
rho1[2][x_, y_, 
          F2_] := {{0, 
            Max[F2[x, y][[2]] - F2[
                    x, y][[1]], 
              0]}, {Max[F2[x, 
                    y][[1]] - F2[x, 
                    y][[2]], 0], 
            0}} ;

PC[x_, y_, F1_, F2_] := 
        Join[Flatten[
            Transpose[rho1[1][x, y, F1]] .  
                x - {{x[[1]], 0}, {0, 
                    x[[2]]}} . rho1[
                    1][x, y, F1] . {1, 1}], 
          Flatten[Transpose[rho1[2][x, y, F2]] .  
                y - {{y[[1]], 0}, {0, 
                    y[[2]]}} . rho1[
                    2][x, y, F2] . {1, 1}]] ;

RPCharacterization[PC] = Nash ;

Projection Dynamic

closetozero = 10^(-7) ;

Projection[x_, y_, F1_, 
          F2_] := {Evaluate[
            If[x[[1]] > 
                  closetozero  ||  (F1[x, 
                        y][[1]] ≥ F1[x, y][[2]]), 
              If[x[[2]] > 
                    closetozero  ||  (F1[x, 
                          y][[2]] ≥  F1[x, 
                          y][[1]]), F1[x, y][[1]] - ((1/2) ((∑_(j = 1)^2) F1[x, y][[
                            j]])), 0], 0]], 
          Evaluate[
            If[x[[2]] > 
                  closetozero  ||  (F1[x, 
                        y][[2]] ≥  F1[x, y][[1]]), 
              If[x[[1]] > 
                    closetozero  ||  (F1[x, 
                          y][[1]] ≥  F1[x, 
                          y][[2]]), F1[x, y][[2]] - ((1/2) ((∑_(j = 1)^2) F1[x, y][[
                            j]])), 0], 0]], 
          Evaluate[
            If[(y[[1]] > 
                    closetozero  ||  (F2[x, 
                          y][[1]] ≥ F2[x, 
                          y][[2]])) , If[(y[[2]] > 
                      closetozero  ||  (F2[x, 
                            y][[2]] ≥  F2[x, 
                            y][[1]])), F2[x, 
                      y][[1]] - ((1/2) ((∑_(j = 1)^2) F2[x, y][[
                            j]])), 0], 0]], 
          Evaluate[
            If[(y[[2]] > 
                    closetozero  ||  (F2[x, 
                          
                          y][[2]] ≥  F2[x, 
                          y][[1]])), If[(y[[1]] > 
                      closetozero  ||  (F2[x, 
                            y][[1]] ≥ F2[x, 
                            y][[2]])) , F2[x, 
                      y][[2]] - ((1/2) ((∑_(j = 1)^2) F2[x, 
                              y][[j]])), 0], 0]] } ;

RPCharacterization[Projection] = Nash ;

Other

Define any other dynamic you want to use. When defining be sure you follow same style with previously defined dynamics. Otherwise, you may need to do major changes in the program.

Other[x_, y_, F1_, 
          F2_] := {x[[1]]  Fhat[
                  1][x, y, F1][[1]], 
          x[[2]]  Fhat[1][x, y, 
                F1][[2]], 
          y[[1]]  Fhat[2][x, y, 
                F2][[1]], 
          y[[2]]  Fhat[2][x, y, 
                F2][[2]]};

RPCharacterization[Other] = Automatic  ;

Choice of contour function

Choose a contour function here.  To use one of the built-in contour functions,  select the rhs of "SimplexContourFunction = xxx" (excluding the semicolon) and then press the button corresponding to the contour function you desire.  If you choose StableGameLyapunov[Dynamic], the program will use the Lyapunov function that is appropriate for the dynamic you selected in the previous section.

({{Speed , L1Speed , NormalFormPotential , StableGameLyapunov[Dynamic] , Other}})

SimplexContourFunction = Speed  ;

Definitions of contour functions

Speed functions

Speed is the standard (Euclidean) speed on the simplex.  L1Speed is self-explanatory.   RSSpeed is the speed of the replicator dynamic after solution trajectories have been moved to the sphere.

Speed[x_, y_] := 
        If[Dynamic ===  Projection, 
          Sqrt[(∑_(i = 1)^4) PhiF[x, y, F1, F2][[i]]^2] , 
          Sqrt[(∑_(i = 1)^4) Dynamic[x, y, F1, F2][[i]]^2] ];

L1Speed[x_, 
          y_] := (∑_(i = 1)^4) 
            Abs[Dynamic[x, y, F1, F2][[
              i]]] ;

Potential functions

NormalFormPotential[x_, y_] := x . (A . y ) ;

Lyapunov functions

NEQ = 
        If[findnashequilibria == 1, NashEq[F1, F2]] ;

StableGameLyapunov[Replicator][x_, 
          y_] := - ((∏_(i = 1)^2) x[[i]]^NEQ[[1, 1, i]]) ((∏_(j = 1)^2) y[[j]]^NEQ[[1, 2, j]])   ;

StableGameLyapunov[Logit[η]][x_, 
          y_] :=  η Log[(∑_(i = 1)^2) (E^(F1[x, y][[i]]/η))] - (x . 
                F1[x, y] - η (∑_(i = 1)^2) 
                      x[[i]] 
                    Log[Max[x[[i]], 
                        10^(-7)]])   +  η Log[(∑_(i = 1)^2) (E^(F2[x, y][[i]]/η))] - (y . 
                F2[x, y] - η (∑_(i = 1)^2) 
                      y[[i]] 
                    Log[Max[y[[i]], 
                        10^(-7)]]) ;

StableGameLyapunov[Projection][x_, 
          y_] :=   (∑_(i = 1)^2)(x[[i]] - 
                  NEQ[[1, 1, 
                    i]])^2   + (∑_(i = 1)^2)(y[[i]] - 
                  NEQ[[1, 2, 
                    i]])^2;

StableGameLyapunov[BNN][x_, 
          y_] := (∑_(i = 1)^2) Fhatplus[1][x, y, F1][[i]]^2   + (∑_(i = 1)^2) Fhatplus[2][x, y, F2][[i]]^2  ;

StableGameLyapunov[PD][x_, 
          y_] := (∑_(i = 1)^2) x[[
                i]] (∑_(j = 1)^2)(((Max[F1[x, y][[j]] - F1[x, y][[i]], 0]))^2)   + (∑_(i = 1)^2) y[[
                i]] (∑_(j = 1)^2)(((Max[F2[x, y][[j]] - F2[x, y][[i]], 0]))^2)   ;

Other contour functions

Define your own contour function here.

Specification of graphical output

Shared parameters for contour function output

These parameters are used for contour function output in all of the types of diagrams described below.  If you do not include contour functions in your output, these parameters are not used.
color = 1 generates color output; color = 0 generates black and white.  The scales used are drawn below.
colorbar = 1 to create a color scale (alternatively, a black and white scale) along with your contour plot.
plotprecision specifies the number of points sampled in each dimension when drawing the contour plots.  Typicaly, this option is the main determinant of the program's running time.  Settings between 50 and 200 are accurate enough for most purposes.  For fast drafts a setting of 20 is adequate.
numberofcontours specifies the number of contour levels shown in the contour plots.
showcontourformula = 1 displays the formula for the contour function in the output notebook.

color = 1 ; 
colorbar = 0 ;
plotprecision = 70 ;
numberofcontours = 50 ;
showcontourformula = 1 ;

Vector field diagram

vectorfielddiagram = 1 draws a vector field diagram on the simplex.  The vector fields should take values in the tangent space of the simplex (i.e., values whose components sum to zero).  
drawoncontour = 1 draws the vector field on top of a contour plot.

vectorfielddiagram = 1 ;
drawoncontour = 0 ;

Instructions for specification of vector fields

Either one or two vector fields is drawn, according to whether drawvectorfield1, drawvectorfield2, or both are set to 1.
vectorfield? specifies the vector field to be drawn.  (Here and below, "?" equals 1 or 2.)  Common specifications for this are Dynamic (the dynamic specified in the "Choice of dynamic" section) and PhiF (the projection of the payoff vector field onto the tangent space of the simplex - see "Payoff-related definitions" in the "Choice of game"section).
scale? rescales the vectors, as is often needed to prevent the vector field arrows from being too long.  The program calibrates the output so that setting scale1 = scale2 makes the longest vector in each vector field have the same length.
vf?thickness is the thickness of the "shafts" of the arrows.
arrowheadlengtsvf? is self-explanatory; set it to zero if you do not want arrowheads.
colorvf? determines the color of the vector field.  It should be set to Hue[t] or GrayLevel[t] using the charts below.

The base points for the vectors in the diagram can be chosen from a grid 
(i) with the grid including points on the boundary (i.e., the dots in the diagram below);
(ii) with the grid excluding points on the boundary (i.e.,the "x"s in the diagram below).
Option (i) is the default setting.
To switch options, select the rhs of "datavf = ..." (excluding the semicolon) and press the appropriate palette button.
gridpoints1dvf specifies the number of initial conditions on each boundary of the grid.
drawrestpoints = 1 causes the rest points of the dynamic to be drawn on the diagram.
arrowsfromrestpoints = 1 causes vector field arrows to be drawn from the rest points (in addition to all other arrows).

drawvectorfield1 = 1 ;
vectorfield1 = PhiF ;
vf1thickness =  .001 ;
scale1 =  .06 ;
arrowheadlengthvf1 =  .01 ;
colorvf1 =  Hue[ .64] ;

drawvectorfield2 = 1 ;
vectorfield2 = Dynamic ;
vf2thickness =  .001 ;
scale2 =  .06 ;
arrowheadlengthvf2 =  .01 ;
colorvf2 =  Hue[ .84] ;

    
gridpoints1dvf =  25 ;

    
drawrestpoints = 1;
arrowsfromrestpoints =  1 ;

({{Grid initial conditions (boundary included)}, {Grid initial conditions (boundary excluded)}})

datavf = 
        Flatten[Table[{{(i/gridpoints1dvf), 
                1 - (( i)/gridpoints1dvf)}, {(j/gridpoints1dvf), 
                1 - (j/gridpoints1dvf)}}, {i, 0, gridpoints1dvf, 1}, {j, 0, 
              gridpoints1dvf, 1}] , 1] ;

Phase diagram

phasediagram = 1 creates a phase diagram on the simplex.  The solution trajectories to be included are specified below.
includecontourplot = 1 draws the phase diagram on top of a contour plot.
ppoints determines how many points will be used when drawing the solution trajectories.

phasediagram = 1 ;
includecontourplot = 1 ;
ppoints = 50 ;

Instructions for "automated" specification of solution trajectories

This section specifies the solution trajectories to be included in the phase diagram using built-in specifications.  These are of three types:
(i) initial conditions drawn from a grid, with the grid including points on the boundary ;  
(ii) initial conditions drawn from a grid, with the grid excluding points on the boundary ;
(iii) randomly chosen initial conditions.  
To use one of these options, select the rhs of "data = ..." (excluding the semicolon) and then press the appropriate palette button.
endtime determines the time interval [0, endtime] of each solution.  
For the grid options, gridpoints1d specifies the number of initial conditions on each boundary of the grid.  
For the random option, randompoints specifies the number of random initial conditions.
numberofarrows specifies how many arrows to draw on each solution.  
arrowheadlength is self-explanatory; it can be set to zero if arrowheads are not desired.  
If numberofarrows = 1 is specified, the arrow appears at the end of each solution trajectory.  If numberofarrows = n, arrows are drawn at times {endtime/n,  2 endtime/n, ... , endtime}.  For more on how arrows are drawn, see "Adjusting positions of arrow heads and 'tails'" below.

Instructions for "manual" specification of solution trajectories

The user can also specify which solution trajectores to graph directly.  To do this, data must be defined as a table of the following form: {{initial condition 1, T1, {color1, thikness1, dashing1}, {Arrow positions 1}, {Head lengths 1}}, {initial condition 2, T2, {color2, thickness2, dashing2}, {Arrow positions 2}, {Head lengths 2}}, ... }
Each initial condition is entered as a list {x2, y2} containing the weight placed on each population's second strategy.  Ti determines time interval [0, Ti] of the ith solution.  Arrow positions i and Head lengths i should be vectors of the same length, where this length is the number of arrowheads to be drawn on the solution trajectory.  Arrow positions are stated as points in the time interval [0, Ti].  For no arrows, set Arrow positions i and Head lengths i to {0}.

{colori, thicknessi, dashingi} specifies the appearance of solutions trajectories.  You can specify any collection of graphics primitives you like (even ones beyond the three mentioned here), and you canplace them in any order.  To use the default appearance, write {}. 
colori can be Hue[?], GrayLevel[?] or any other color command;
thicknessi can be Thickness[?] or AbsoluteThickness[?];
dashingi is specified as Dashing[{???}].

endtime = 1 ; 

    gridpoints1d = 9; 

    randompoints = 60 ;

    

    numberofarrows = 2 ;

    arrowheadlength =  .015 ;

({{Grid initial conditions (boundary included)}, {Grid initial conditions (boundary excluded)}, {Random initial conditions}})

data = 
        Flatten[Table[{{1 - 
                  (1/(2 gridpoints1d)) - (( i)/gridpoints1d), 
                1 - (1/(2 gridpoints1d)) - (( j)/gridpoints1d)}, 
              endtime, {}, 
              Table[(((numberofarrows - i + 1) endtime)/numberofarrows), {i, numberofarrows}], 
              Table[arrowheadlength, {i, numberofarrows}]}, {i, 0, 
              gridpoints1d - 1, 1}, {j, 0, gridpoints1d - 1, 1}] , 
          1]   ;

Tuning rest point search

RestSearchInit contains the initial conditions used in the FindRoot command when the program searches for rest points.  If you suspect that the program has missed some rest points, add initial conditions in the regions where the missing rest points should be.  Each initial condition is written as {y2, x1}; x2 and y1 are computed as 1 - x1 and 1 -y2.
RestTolerance is a parameter in a filtering procedure used to eliminate false rest points.  If the speed of the dynamic at a prospective rest point is above (RestTolerance) * (the maximum speed of the dynamic), the prospective rest point is thrown out.

RestSearchInit = {{0.001, 0.001}, {0.001, 0.101}, {0.001, 
            0.201}, {0.001, 0.301}, {0.001, 0.401}, {0.001, 0.501}, {0.001, 
            0.601}, {0.001, 0.701}, {0.001, 0.801}, {0.001, 0.901}, {0.001, 
            0.999}, {0.101, 0.001}, {0.101, 0.101}, {0.101, 0.201}, {0.101, 
            0.301}, {0.101, 0.401}, {0.101, 0.501}, {0.101, 0.601}, {0.101, 
            0.701}, {0.101, 0.801}, {0.101, 0.901}, {0.099, 0.999}, {0.201, 
            0.001}, {0.201, 0.101}, {0.201, 0.201}, {0.201, 0.301}, {0.201, 
            0.401}, {0.201, 0.501}, {0.201, 0.601}, {0.201, 0.701}, {0.201, 
            0.801}, {0.201, 0.901}, {0.199, 0.999}, {0.301, 0.001}, {0.301, 
            0.101}, {0.301, 0.201}, {0.301, 0.301}, {0.301, 0.401}, {0.301, 
            0.501}, {0.301, 0.601}, {0.301, 0.701}, {0.301, 0.801}, {0.301, 
            0.901}, {0.299, 0.999}, {0.401, 0.001}, {0.401, 0.101}, {0.401, 
            0.201}, {0.401, 0.301}, {0.401, 0.401}, {0.401, 0.501}, {0.401, 
            0.601}, {0.401, 0.701}, {0.401, 0.801}, {0.401, 0.901}, {0.399, 
            0.999}, {0.501, 0.001}, {0.501, 0.101}, {0.501, 0.201}, {0.501, 
            0.301}, {0.501, 0.401}, {0.501, 0.501}, {0.501, 0.601}, {0.501, 
            0.701}, {0.501, 0.801}, {0.501, 0.901}, {0.499, 0.999}, {0.601, 
            0.001}, {0.601, 0.101}, {0.601, 0.201}, {0.601, 0.301}, {0.601, 
            0.401}, {0.601, 0.501}, {0.601, 0.601}, {0.601, 0.701}, {0.601, 
            0.801}, {0.601, 0.901}, {0.599, 0.999}, {0.701, 0.001}, {0.701, 
            0.101}, {0.701, 0.201}, {0.701, 0.301}, {0.701, 0.401}, {0.701, 
            0.501}, {0.701, 0.601}, {0.701, 0.701}, {0.701, 0.801}, {0.701, 
            0.901}, {0.699, 0.999}, {0.801, 0.001}, {0.801, 0.101}, {0.801, 
            0.201}, {0.801, 0.301}, {0.801, 0.401}, {0.801, 0.501}, {0.801, 
            0.601}, {0.801, 0.701}, {0.801, 0.801}, {0.801, 0.901}, {0.799, 
            0.999}, {0.901, 0.001}, {0.901, 0.101}, {0.901, 0.201}, {0.901, 
            0.301}, {0.901, 0.401}, {0.901, 0.501}, {0.901, 0.601}, {0.901, 
            0.701}, {0.901, 0.801}, {0.901, 0.901}, {0.899, 0.999}, {0.999, 
            0.001}, {0.999, 0.099}, {0.999, 0.199}, {0.999, 0.299}, {0.999, 
            0.399}, {0.999, 0.499}, {0.999, 0.599}, {0.999, 0.699}, {0.999, 
            0.799}, {0.999, 0.899}, {0.999, 0.999}};

RestTolerance =  .01 ;

Tuning stability tests

The first four parameters determine how rest points are checked for stability.  The parameter numtest specifies how many trajectories will be checked.  These initial conditions are random points at distance eps from the rest point.  Solutions from these points over the time interval [0, bigT] are computed.  If any of these solutions winds up more than delta away from the rest point, the rest point is considered unstable; otherwise it is considered stable.

numtest = 5 ; 

    eps =  .01 ; 

    delta =  .01 ;

    bigT = 1 ;

Adjusting positions of arrow heads and "tails"

To draw an arrow, the program must know where to put the tail of the arrow and where to put the head.  Even if one does not want the base of the arrow to appear in the picture, its position of the tail must be specified so that the arrow points in the right direction.  
The parameter relarrowbase below determines the positions of arrows' tails relative to their heads when arrows are drawn on solution trajectories.  In particular, if an arrow is supposed to appear at a solution trajectory's position at time t, the program draws the head at this point, placing the tail at the solution trajectory's position at time (t – endtime/relarrowbase).  
A value of .008 for relarrowbase works well in most cases:  the arrow base is not visible, and the arrow points in the right direction.  If arrow bases are visible, a smaller value of relarrowbase should be specified.  If arrows point in the wrong directions, a larger value of relarrowbase should be chosen.
The parameter arrowheadcenter  specifies the location of the center of the arrowheads' bases as a fraction of the arrowheads' lengths.

relarrowbase =  .005 ;
arrowheadcenter =  .8 ;

Adding extra graphic elements

This section allows the user to add dots, lines, arrows, text, and color solution trajectories to a phase diagram/contour plot.

Enter the data for dots as follows:  dotdata = {{center1, radius1, color1, open1}, {center2, radius2, color2, open2}, ...}. color1 should be a Mathematica-defined color function (Hue[.], GrayLevel[.], etc).  Set open1 = 0 for solid dot and = 1 for a dot with a white center. centers should be entered like {{x1,x2},{y1,y2}}.
 
Enter the data for lines as follows:  linedata = {{startingpoint1, endpoint1, thickness1, dashing1, color1, arrowheadlength1}, {startingpoint2, endpoint2, thickness2, dashing2, color2, arrowheadlength2},...}.  If you want an arrow at the end of the line enter the value for arrowheadlength specifying the size of arrow head; if you do not want an arrow simply set it to 0. Thickness and Dashing are usual Mathematica commands. starting points should be entered like {{x1,x2},{y1,y2}}. endpoints should be entered like {{x1,x2},{y1,y2}}. Line start at startingpoint and ends at the endpoints as names suggest.

Enter data for text as follows:  textdata = {{"label1", position1, "fontfamily1","fontslant1", fontsize1}, {"label2", position2, "fontfamily2","fontslant2", fontsize2}, ...}. Remember that you are restricted to positions in PlotRange -> {{-.08, 1.1}, {-.08, 1.1}}.  positions should be entered in rectangular coordinates {x,y}.

Enter data for extra solution trajectories as follows:  cycledata = {{initial condition on simplex 1, endtime1 for NDSolve, beginningtime1 to be used in in the plot, {color1, thickness1}, {arrow position1}, {head lengths1}}, {initial condition on simplex 2, endtime2 for NDSolve, beginningtime2 to be used in in the plot, {color2, thickness2}, {arrow position2}, {head lengths2}},...}. 
The program computes a solution to the dynamic from time 0 through endtime?, and then plots the solution trajectory from the beginingtime? through endtime.   color? should be a Mathematica-defined color function (Hue[.], GrayLevel[.], etc).  (both color and thickness can be ignored, just write {} if you do not want to bother with these.  arrowpositon? must be a number strictly greater than beginingtime? and less than endtime (it must be strictly bigger than begining time).  Each trajectory can have multiple arrows, but even if there is only one arrow, be sure to put brackets around the arrow positions and head lengths. Initial conditions should be entered like {x2,y2}.

Write dotdata = {},  linedata = {}, textdata = {}, and cycledata = {} if you do not want to add anything.  (This is the default setting.)

Example:  This example generates one line, one dot, and one piece of text, one special solution trajectory (cycledata).
dotdata={{{{1/2,1/2},{1/2,1/2}}, .01, GrayLevel[0], 1}} ;    [draws a black dot with white center at the center of the simplex]
linedata={{{{1/2, 1/2},{1,0}},{{1/2, 1/2},{0,1}}, Thickness[.004], Dashing[{}], GrayLevel[0], 0}} ;   [draws a horizantal line through at y=1/2]
textdata={{"trial", {1/2,1/2}, "Palatino-Roman", "Plain", 12}} ; 
cycledata={{{8,.9}, 50, 45, {Hue[.83], Thickness[.005]}, {47}, {.03}}} ;

dotdata = {} ;
linedata = {} ;
textdata =  {} ;
cycledata = {} ;

uddots = 
        Graphics[
          Table[{dotdata[[i, 3]], 
              Disk[{dotdata[[i, 1, 2, 
                    2]], 
                  dotdata[[i, 1, 1, 
                    1]]}, 
                dotdata[[i, 
                  2]]]}, {i, Dimensions[
                  dotdata][[1]]}]] ;
udwdots = 
        Graphics[
          DeleteCases[
            Table[If[
                dotdata[[i, 4]] == 
                  1, {GrayLevel[1], 
                  Disk[{dotdata[[i, 1, 2, 
                        2]], 
                      dotdata[[i, 1, 1, 
                        1]]},  .6 dotdata[[i, 2]]]}], {i, Dimensions[
                    dotdata][[1]]}], 
            Null]] ;
extralines = 
        Graphics[
          DeleteCases[
            Table[{linedata[[i, 3]], 
                linedata[[i, 4]], linedata[[i, 5]], 
                Arrow[{linedata[[i, 1, 2, 
                      2]], 
                    linedata[[i, 1, 1, 
                      1]]}, {linedata[[
                      i, 2, 2, 2]], 
                    linedata[[i, 2, 1, 
                      1]]}, HeadCenter ->  4/5, 
                  HeadLength -> 
                    linedata[[i, 
                      6]]]}, {i, Dimensions[
                    linedata][[1]]}], 
            Null]] ;
extratext = 
        Graphics[
          DeleteCases[
            Table[{Text[
                  textdata[[i, 1]], 
                  textdata[[i, 2]], 
                  TextStyle -> {FontFamily -> 
                        textdata[[i, 
                          3]], 
                      FontSlant -> 
                        textdata[[i, 
                          4]], 
                      FontSize -> 
                        textdata[[i, 
                          5]]}]}, {i, Dimensions[
                    textdata][[1]]}], Null]];

3D graph of the contour function

graphcontour3d = 1 creates a 3d plot of the contour function. 
shiftpercentage specifies the distance above the simplex that the lowest point on the contour function will appear, as a percentage of the distance from the minimum value of the function to the maximum value of the function.
grayline = 1 draws a gray line on the edges of the 3d plot.
heightratio specifies the ratios of the three coordinates in the graphical output.  The default setting is {1, 1, 1.25}.
  
vpoint determines the value for the Mathematica command ViewPoint and hence the point from which the 3d plot is observed.  In Mathematica, ViewPoint -> {x, y, z} gives the position of the view point relative to the center of the three-dimensional box that contains the object being plotted.  The view point is given in a special scaled coordinate system in which the longest side of the bounding box has length 1.  (When heightratio is at its default setting of {1, 1, 1.25}, the longest side corresponds to the z coordinate.  The center of the bounding box is taken to have coordinates {0, 0, 0}. Choosing ViewPoint further away from the object reduces the distortion associated with perspective.

graphcontour3d = 0 ;

    
grayline = 1 ;
shiftpercentage =  .2 ;
heightratio = {1, 1, 1.25} ;
vpoint = {-2.730,  -1.673,  1.096} ;

Program Core

Solving the differential equation

time1 = TimeUsed[] ;

x[t] = {x1[t], 1 - x1[t]} ;
y[t] = {y1[t], 1 - y1[t]} ;

G=Table[NDSolve[{x1'[t]==Dynamic[x[t], y[t], F1, 
                          F2][[1]],y1'[t]==Dynamic[x[t], y[t], F1, 
                          F2][[3]],x1[0] == 
                      1 - data[[i, 
                            1]][[1]],y1[0] == 
                      1 - data[[i, 
                            1]][[2]]},{x1, y1},{t, 0, 
                  data[[i, 2]]}],{i, Dimensions[
                  data][[1]]}] ;

GCycle=If[Dimensions[
                  cycledata][[1]] == 0,{},Table[NDSolve[{x1'[t]==Dynamic[x[t], y[t], F1, 
                              F2][[1]],y1'[t]==Dynamic[x[t], y[t], F1, 
                              F2][[3]],x1[0] == 
                          1 - cycledata[[i, 
                                1]][[1]],y1[0] == 
                          1 - cycledata[[i, 
                                1]][[2]]},{x1, y1},{t, 0, 
                      cycledata[[i, 
                        2]]}],{i, Dimensions[
                      cycledata][[1]]}]] ;

Plotting the solutions

If[vectorfielddiagram == 0,  drawvectorfield1 = 0 ;  
        drawvectorfield2 = 0 ] ;

    
If[phasediagram == 0, includecontourplot = 0] ;

Rest points and stability test

TolLevel = 
        Times[RestTolerance,  
          Max[Flatten[
              Table[Chop[Speed[{i, 1 - i}, {j, 1 - j}], (10^(-5))], {i, 0, 
                  1,  .05}, {j, 0, 1,  .05}]]]]  ;

If[findrestpoints ==  1, 
        RePos = RestPoints[Dynamic, F1, F2, RestSearchInit, TolLevel] ; 
        TestStability[Dynamic, F1, F2, bigT,  numtest,  eps,  delta, 
          RestSearchInit, TolLevel]] ;

Drawing vector fields on the simplex

Needs["Graphics`PlotField`"] ;

If[
        vectorfielddiagram == 
            1  &&  (drawvectorfield1 == 1  ||  
              drawvectorfield2 == 1),  
colvf1[
            dum_] = colorvf1 ; 
colvf2[dum_] = 
          colorvf2 ; 

If[
          drawvectorfield1 == 0, vectorf1 = Graphics[{}], 
          vectorf1 = 
            ListPlotVectorField[ 
              If[findrestpoints == 1  &&  
                  arrowsfromrestpoints == 1, 
                Join[Table[{{datavf[[i, 2, 2]], 
                        datavf[[i, 1, 1]]}, {vectorfield1[
                            datavf[[i, 1]], datavf[[i, 2]], F1, 
                            F2][[4]], vectorfield1[datavf[[i, 1]], datavf[[i, 2]], F1, 
                            F2][[1]]}}, {i, Dimensions[datavf][[1]]}],  
                  Table[{{testrest[[i, 2, 2]], 
                        testrest[[i, 1, 1]]}, {vectorfield1[
                            testrest[[i, 1]], testrest[[i, 2]], F1, 
                            F2][[4]], vectorfield1[testrest[[i, 1]], testrest[[i, 2]], F1, 
                            F2][[1]]}}, {i, Dimensions[testrest][[1]]}] ],  
                Table[{{datavf[[i, 2, 2]], 
                      datavf[[i, 1, 1]]}, {vectorfield1[
                          datavf[[i, 1]], datavf[[i, 2]], F1, 
                          F2][[4]], vectorfield1[datavf[[i, 1]], datavf[[i, 2]], F1, 
                          F2][[1]]}}, {i, Dimensions[datavf][[1]]}] ], ScaleFactor ->  scale1 , 
              HeadCenter -> arrowheadcenter, 
              HeadLength ->  arrowheadlengthvf1,  
              ColorFunction ->  ({colvf1[#] , 
                      Thickness[vf1thickness]} &),  
              DisplayFunction ->  
                Identity ] ] ; 

If[
          drawvectorfield2 == 0, 
          vectorf2 = Graphics[{}], 
vectorf2 = 
            ListPlotVectorField[ 
              If[findrestpoints == 1  &&  
                  arrowsfromrestpoints == 1, 
                Join[Table[{{datavf[[i, 2, 2]], 
                        datavf[[i, 1, 1]]}, {vectorfield2[
                            datavf[[i, 1]], datavf[[i, 2]], F1, 
                            F2][[4]], vectorfield2[datavf[[i, 1]], datavf[[i, 2]], F1, 
                            F2][[1]]}}, {i, Dimensions[datavf][[1]]}],  
                  Table[{{testrest[[i, 2, 2]], 
                        testrest[[i, 1, 1]]}, {vectorfield2[
                            testrest[[i, 1]], testrest[[i, 2]], F1, 
                            F2][[4]], vectorfield2[testrest[[i, 1]], testrest[[i, 2]], F1, 
                            F2][[1]]}}, {i, Dimensions[testrest][[1]]}] ],  
                Table[{{datavf[[i, 2, 2]], 
                      datavf[[i, 1, 1]]}, {vectorfield2[
                          datavf[[i, 1]], datavf[[i, 2]], F1, 
                          F2][[4]], vectorfield2[datavf[[i, 1]], datavf[[i, 2]], F1, 
                          F2][[1]]}}, {i, Dimensions[datavf][[1]]}] ], ScaleFactor ->  scale2 , 
              HeadCenter -> arrowheadcenter, 
              HeadLength ->  arrowheadlengthvf2,  
              ColorFunction ->  ({colvf2[#] , 
                      Thickness[vf2thickness]} &),  
              DisplayFunction ->  Identity ] ] ] ;

Drawing Solutions to Diff Eqs if there is no contour plot

e = 
        Table[ParametricPlot[
            Evaluate[{1 - y1[t], 
                  x1[t]} /. G[[
                  i]] ], {t, 0, data[[i, 2]]},  
            DisplayFunction  ->  Identity, 
            PlotPoints ->  ppoints,  Axes -> {False, False}, 
            PlotRange -> {{- .08, 1.08}, {- .08, 1.08}}, 
            PlotStyle ->  {data[[i, 
                  3]]}, 
            AspectRatio -> 
              Automatic], {i, Dimensions[
                data][[1]]}] ;

ecycle = 
        If[Dimensions[
                cycledata][[1]] == 0, Graphics[{}], 
          Table[ParametricPlot[
              Evaluate[{1 - y1[t], x1[t]} /. 
                  GCycle[[i]]], {t, 
                cycledata[[i, 3]], 
                cycledata[[i, 2]]}, 
              DisplayFunction -> Identity, Axes -> {False, False}, 
              PlotRange -> {{- .08, 1.08}, {- .08, 1.08}}, 
              PlotStyle -> {cycledata[[i, 
                    4]]}, 
              AspectRatio -> 
                Automatic], {i, Dimensions[
                  cycledata][[1]]}]];

Specifying color scale/gray scale

The approximate minimum and maximum values of the contour function are computed here.  These values are then used to assign colors/gray levels optimally.

If[
        includecontourplot == 1  || drawoncontour == 1 || 
          graphcontour3d  ==  1 ,  
        If[SimplexContourFunction  ===  Speed, MinAttained = 0, 
          nrmlzdmin = 
            Flatten[Table[
                Chop[SimplexContourFunction[{i, 1 - i}, {j, 1 - j}], 
                  (10^(-5))], {i, 0, 1,  .025}, {j, 0, 
                  1,  .025}]] ; 
MinAttained = 
            Min[nrmlzdmin]]] ;

    
If[
        includecontourplot == 1  ||  drawoncontour == 1 || 
          graphcontour3d  ==  1, 
        nrmlzd = 
          Flatten[Table[
              Chop[SimplexContourFunction[{i, 1 - i}, {j, 1 - j}], 
                (10^(-5))], {i, 0, 1,  .025}, {j, 0, 
                1,  .025}]] ; 
MaxAttained = 
          Max[nrmlzd] ] ;

    
distancefrommin = 
        shiftpercentage (MaxAttained - 
              MinAttained) ;

    
btwf[Ψ_] = 
        1 - Max[0, 
            Min[((Ψ - MinAttained)/(MaxAttained - MinAttained)), 1]] ;

    
If[(includecontourplot == 1  || 
              drawoncontour == 1 ||    graphcontour3d  ==  
                1)  &&  color == 0, 
        ConColFun[Ψ_] = 
          GrayLevel[ .(95/1.45) (1 - 
                  btwf[Ψ] +  .45) ]] ;
If[(includecontourplot == 1  || 
              drawoncontour == 1 ||  
              graphcontour3d  ==  1) &&  color == 1, 
        ConColFun[Ψ_] = 
          Hue[ Which[
              btwf[Ψ] ≤ (1/2), (2/5) 
                btwf[Ψ], 
              (1/2) < btwf[Ψ] ≤ 
                (2/3),  (6/5) btwf[Ψ] - (2/5), 
              (2/3) < btwf[Ψ], (3/4) btwf[Ψ] - 
                (1/10)]]] ;

    
If[(includecontourplot == 1  ||  
              drawoncontour == 1) &&  colorbar == 1, 
        simplexcolorbar = 
          If[color == 
              1, {Graphics[({Hue[
                          Which[1 - #1  ≤ 
                              (1/2), (2/5) (1 - #1), 
                            (1/2) < (1 - #1) ≤ 
                              (2/3),  (6/5) (1 - #1) - (2/5), 
                            (2/3) < (1 - #1), (3/4) (1 - #1) - 
                              (1/10)]], 
                        Rectangle[{1.03, ((#1 + 1.5)/3)}, {1.08, ((#1  + 1.51)/3)}]} &) /@ Range[0, 1,  .01]],  
              Graphics[{Text[
                    PaddedForm[ MinAttained,  {3, 2}], {1.08,  .48}],  
                  Text[PaddedForm[
                      MaxAttained, {3, 
                        2}], {1.08,  .85}]}]}, {Graphics[({GrayLevel[ .(95/1.45) ((#1) +  .45)], 
                        Rectangle[{1.03, ((#1 + 1.5)/3)}, {1.08, ((#1  + 1.51)/3)}]} &) /@ Range[0, 1,  .01]],  
              Graphics[{Text[
                    PaddedForm[MinAttained, {3, 2}], {1.08,  .48}],  
                  Text[PaddedForm[ 
                      MaxAttained, {3, 2}], {1.08,  .85}]}]}] , 
        simplexcolorbar = Graphics[{}] ] ;

    
If[includecontourplot == 1, 
        simplexcolorbarph = simplexcolorbar, 
        simplexcolorbarph = Graphics[{}]] ;
If[drawoncontour == 1, simplexcolorbarvf = simplexcolorbar, 
        simplexcolorbarvf = Graphics[{}]] ;

Drawing Solutions to Diff Eqs if there is a contour plot

If[
        includecontourplot == 1  ||  
          drawoncontour == 1, 
     g[y_] := 
          SimplexContourFunction[{y[[2]], 1 - y[[2]]}, {1 - 
                y[[1]], 
              y[[1]]}] ;  

con = 
          ContourPlot[g[{y1, y2}], {y1, 0, 1}, {y2, 0, 1}, 
            DisplayFunction  ->  Identity,  
            PlotPoints -> plotprecision, 
            Contours -> numberofcontours, Frame -> None, 
            AspectRatio -> Automatic,  ContourLines  ->  False, 
            ColorFunctionScaling  ->  False ,  
            ColorFunction ->  (ConColFun[#] &)] , 
        con = Graphics[{}]] ;

    
If[includecontourplot == 1 , conph = con, 
        conph = Graphics[{}]] ;
If[drawoncontour == 1, convf = con, 
        convf = Graphics[{}]] ;

Drawing 3D Plot of contour function

\(\(If[
        graphcontour3d == 1, \ 
frestricted3d[x_, y_] := 
          If[x[\([1]\)] ≥ 0 && x[\([2]\)] ≥ 0 && 
              y[\([1]\)] ≥ 0 && y[\([2]\)] ≥ 0 && 
              x[\([1]\)] + x[\([2]\)] == 1 && 
              y[\([1]\)] + y[\([2]\)] == 1, 
            SimplexContourFunction[x, y], \(-\ 
              Infty\)]\ ; 


        g3d[y_] := 
          frestricted3d[{y[[2]], 
              1 - y[[2]]}, {1 - 
                y[[1]], 
              y[[1]]}]\ ; \ 

confunplt = \ 
          Plot3D[\ g3d[{y1, y2}]\ , \ {y1, 0, 1}, {y2, 0, 1}, 
            DisplayFunction\  -> \ 
              Identity, 
PlotPoints\  -> \ 
              plotprecision, \ 
            PlotRange\  -> \ {MinAttained - distancefrommin, 
                MaxAttained +  .05}, \ ViewPoint -> vpoint, \ 
            Lighting -> \ False, 
            ColorFunctionScaling -> \ False\ , \ 
            ColorFunction -> \ \((ConColFun[#] &)\), \ 
Axes\  -> \ False, \ Boxed\  -> \ False, \ 
            Mesh\  -> \ False, \ 
            BoxRatios\  -> \ {1, 1, 
                1.25}]\ ; 
edgeline1 = \ 
          ParametricPlot3D[\ {y1, 0, 
              SimplexContourFunction[{0, 1}, {1 - y1, y1}], 
              GrayLevel[ .3]}\ , \ {y1, 0, 1}, 
            DisplayFunction\  -> \ 
              Identity, 
PlotPoints\  -> 200, \ 
            PlotRange\  -> \ {MinAttained - distancefrommin, 
                MaxAttained +  .05}, \ 
            ViewPoint -> vpoint, \ 
Axes\  -> \ 
              False, \ Boxed\  -> \ 
              False, \ \ BoxRatios\  -> \ {1, 1, 1.25}]\ ; 
        edgeline2 = \ 
          ParametricPlot3D[\ {y1, 1, 
              SimplexContourFunction[{1, 0}, {1 - y1, y1}], 
              GrayLevel[ .3]}\ , \ {y1, 0, 1}, 
            DisplayFunction\  -> \ 
              Identity, 
PlotPoints\  -> \ 200, \ 
            PlotRange\  -> \ {MinAttained - distancefrommin, 
                MaxAttained +  .05}, \ 
            ViewPoint -> vpoint, \ 
Axes\  -> \ 
              False, \ Boxed\  -> \ 
              False, \ \ BoxRatios\  -> \ {1, 1, 
                1.25}]\ ; \ \ edgeline3 = \ 
          ParametricPlot3D[\ {0, y2, g3d[{0, y2}], GrayLevel[ .3]}\ , \ {y2, 
              0, 1}, DisplayFunction\  -> \ 
              Identity, 
PlotPoints\  -> \ 200, \ 
            PlotRange\  -> \ {MinAttained - distancefrommin, 
                MaxAttained +  .05}, \ 
            ViewPoint -> vpoint, \ 
Axes\  -> \ 
              False, \ Boxed\  -> \ False, \ \ BoxRatios\  -> \ 
              heightratio]; \ 
edgeline4 = \ 
          ParametricPlot3D[\ {1, y2, g3d[{1, y2}], GrayLevel[ .3]}\ , \ {y2, 
              0, 1}, DisplayFunction\  -> \ 
              Identity, 
PlotPoints\  -> \ 200, \ 
            PlotRange\  -> \ {MinAttained - distancefrommin, 
                MaxAttained +  .05}, \ 
            ViewPoint -> vpoint, \ 
Axes\  -> \ 
              False, \ Boxed\  -> \ False, \ \ BoxRatios\  -> \ 
              heightratio];\ ];\)\)

Drawing Triangle and Dots

The cells belongs to this group is for drawing white area around the simplex.

Tri = 
        Graphics[
          Line[{{0, 0}, {0, 1}, {1, 1}, {1, 0}, {0, 0}}]];

This part draws the dots for the rest points.  TestStability[Dynamic,F,bigT, numtest, eps, delta,RestSearchInit,TolLevel] must be called in order to determine stability of points.  Hence you may need to modify this part if  StabilityCheck do not work for your choice of dynamic.

If[
        findrestpoints ==  1, 
DOTS = 
          Graphics[
            Table[Disk[{testrest[[i, 2, 
                    2]], 
                  testrest[[i, 1, 1]]}, 
                0.01], {i, Dimensions[
                    testrest][[1]]}]] ; 
WDots = 
          Graphics[
            Table[{GrayLevel[1], 
                Disk[{unstbl[[i, 2, 2]], 
                    unstbl[[i, 1, 1]]}, 
                  0.006]}, {i, Dimensions[
                    unstbl][[1]]}]], 
        DOTS = Graphics[{}] ;  WDots = Graphics[{}] ] ;

These cells draw the triangle, lines, and dots for the 3D plot of contour function

If[graphcontour3d == 1,  
        graphf[x_, y_]  :=  {y[[2]], 
            x[[1]], 
            SimplexContourFunction[x, y]} ; 
graphmin[x_, 
            y_]  :=  {y[[2]], 
            x[[1]], 
            MinAttained - distancefrommin} ; 
Tri3d = 
          Graphics3D[
Line[{{0, 0, 
                  MinAttained - distancefrommin},  {1, 0, 
                  MinAttained - distancefrommin}, {1, 1, 
                  MinAttained - distancefrommin}, {0, 1, 
                  MinAttained - distancefrommin},  {0, 0, 
                  MinAttained - 
                    distancefrommin} } ]] ; 
If[
          findrestpoints == 1, 
          SeveralLines = 
Graphics3D[ {Dashing[{ .01,   .01}],  Flatten[Table[
                    Line[{graphmin[
                          RePos[[i, 1]], 
                          RePos[[i, 
                            2]]], 
                        graphf[
                          RePos[[i, 1]], 
                          RePos[[i, 
                            2]]]}], {i, 
                      1, Dimensions[
                          RePos][[1]], 
                      1}]]}] ; 

SeveralPoints = 
            Graphics3D[ {PointSize [ .02],  
                Flatten[
                  Table[Point[
                      graphmin[
                        stblpnts[[i, 
                          1]], 
                        stblpnts[[i, 
                          2]]]], {i, Dimensions[
                          stblpnts][[1]]}]], GrayLevel[ .5], 
                Flatten[
                  Table[Point[
                      graphmin[
                        unstbl[[i, 1]], 
                        unstbl[[i, 
                          2]]]], {i, Dimensions[
                          unstbl][[1]]}]]}]]] ;

Arrows

ar = 
        Graphics[
          Flatten[Table[{If[Dimensions[
                        data[[i, 
                          3]]][[1]] ≠  
                    0, Join[
                      Select[data[[i, 
                          3]], (MatchQ[#1, Hue[__]]  || 
                              MatchQ[#1, GrayLevel[__]] || 
                              MatchQ[#1, RGBColor[__]] || 
                              MatchQ[#1, CMYKColor[__]]) &], {GrayLevel[
                          0]}][[1]], 
                  GrayLevel[0]], 
                Table[Arrow[
                    Flatten[
                      Evaluate[{1 - 
                              y1[data[[i, 
                                      4]][[j]] - 
                                  relarrowbase data[[i, 
                                      2]]], 
                            x1[data[[i, 
                                    4]][[j]] - 
                                relarrowbase data[[i, 
                                    2]]]} /. G[[i]]]], 
                    Flatten[
                      Evaluate[{1 - 
                              y1[data[[i, 
                                    4]][[j]]], 
                            x1[data[[i, 
                                  4]][[
                                j]]]} /. G[[i]]]], 
                    HeadCenter -> arrowheadcenter, 
                    HeadLength -> data[[i, 
                          5]][[
                        j]]], {j, Dimensions[
                        data[[i, 
                          5]]][[1]]}]}, {i, Dimensions[
                    data][[1]]}]]] ;

arcycle = 
        If[Dimensions[
                cycledata][[1]] == 0, Graphics[{}], 
          Graphics[
            Flatten[Table[{If[Dimensions[
                          cycledata[[i, 
                            4]]][[1]] ≠  
                      0 , Join[
                        Select[
                          cycledata[[i, 
                            4]], (MatchQ[#1, Hue[__]]  || 
                                MatchQ[#1, GrayLevel[__]] || 
                                MatchQ[#1, RGBColor[__]] || 
                                MatchQ[#1, CMYKColor[__]]) &], {GrayLevel[
                            0]}][[1]], 
                    GrayLevel[0]], 
                  Table[Arrow[
                      Flatten[
                        Evaluate[{1 - 
                                y1[cycledata[[i, 
                                        5]][[j]] - 
                                    relarrowbase data[[i, 
                                        2]]], 
                              x1[cycledata[[i, 
                                      5]][[j]] - 
                                  relarrowbase data[[i, 
                                      2]]]} /. 
                            GCycle[[
                              i]]]], 
                      Flatten[
                        Evaluate[{1 - 
                                y1[cycledata[[i, 
                                      5]][[j]]], 
                              x1[cycledata[[i, 
                                    5]][[j]]]} /. 
                            GCycle[[
                              i]]]], 
                      HeadCenter -> arrowheadcenter, 
                      HeadLength -> cycledata[[i, 
                            6]][[
                          j]]], {j, Dimensions[
                          cycledata[[i, 
                            6]]][[1]]}]}, {i, Dimensions[
                      cycledata][[1]]}]]]];

Generating the output notebook

\(\(\( (*\(\(--\(--\(--\(--\(--Redirecting\)\)\)\)\)\ output\ to\ a\ new\ \(\(\(notebook--\)--\)--\)\)\(-\)*) \)\(
\)\(\($nb = 
        NotebookCreate[
          WindowSize -> \ {650, 550}]\ ;\)

    \($DisplayFunction = 
        NotebookWrite[$nb, \ 
            Cell[GraphicsData["PostScript", \ 
                DisplayString[#]], "Graphics"]] &\ ;\)

    \(\($Post = 
        If[# =!= \ Null, 
            NotebookWrite[$nb, 
              Cell[BoxData[
                  MakeBoxes[#, 
                    StandardForm]], "Output"]]] &\ ;\)\(
\)
    \)

    \(Unprotect[Print]\ ;\)

    \(Print[x_]\  := \ 
        NotebookWrite[$nb, 
          Cell[BoxData[MakeBoxes[x, StandardForm]]]]\ ;\)

    \(Print[x__]\  := Map[Print, {x}]\ ;\)

    \(\(Protect[Print]\ ;\)\(
\)
    \)

    \(Unprotect[StylePrint]\ ;\)

    \(StylePrint[x_]\  := \ 
        NotebookWrite[$nb, 
          Cell[BoxData[ToBoxes[x]], "Output"]]\ ;\)

    \(StylePrint[x_, rest___]\  := \ 
        NotebookWrite[$nb, 
          Cell[BoxData[ToBoxes[x]], rest]]\ ;\)

    \(Protect[
        StylePrint];\)

 (*\(--\(--\(--\(--\(--\(--\(--\(--Redirecting\)\)\)\)\)\)\)\)\ code\ \(\(\(\(\(ends--\)--\)--\)--\)--\) - \(--\(---\)\)*) 

    \)\)\)

NotebookWrite[$nb, 
        Cell[BoxData[
            ToBoxes["💡💡💡 DYNAMO \
💡💡💡"]], "Title",  
          FontColor ->  Hue[ .19], Background ->  Hue[ .72], 
          TextAlignment ->  Center, 
          CellFrame ->  True]] ;

Text output - the game

\(\(\(Print[
        StyleForm[
          ColumnForm[{"", "The Game"}], "Section"]]\ ;\)\(
\)
    \)\)
If[nonlineargame  ==  0, 
        Print[SequenceForm["   Normal form game:  ", 
            GridBox[{{""}, {strategy[1]}, {strategy[2]}}]  // 
              DisplayForm , 
GridBox[
                Join[{{strategy[3], strategy[4]}}, U]] // 
              DisplayForm]]] ;

    
Print[
        SequenceForm["   Payoff vector field:  ", "F", {x, y} // 
            MatrixForm,  " = ", 
          Flatten[F[{x1, x2}, {y1, y2}]]  //  
            MatrixForm]] ;

    
If[findnashequilibria == 1,  
        If[  Dimensions[NashEq[F1, F2]][[1]] == 1, 
          Print["   Nash equilibrium:"], 
          Print["   Nash equilibria:"]] ; 
Print[
          SequenceForm["     ", 
            ColumnForm[
              Table[PaddedForm[NashEq[F1, F2][[i]], {7, 7}, 
                  NumberPadding -> {" ", " "}], {i, 
                  1, Dimensions[
                      NashEq[F1, F2]][[1]]}]] ]] ];

Text output - the dynamic

\(\(\(\(\(Print[
        StyleForm[
          ColumnForm[{"", "The Dynamic"}], "Section"]]\ ;\)\(
\)
    \)

    \(\(Print[
        SequenceForm["   Name of dynamic:  ", \ \ Dynamic\ ]]\ ;\)\(
\)
    \)

    \(Print["   Law of motion: ", \ 
        SequenceForm[\ \*""   (Subscript[V, F])"", {x, y} // 
            MatrixForm, \ " = ", 
          Simplify[Dynamic[{x1, x2}, {y1, y2}, F1, F2]] // 
            MatrixForm\ ]\ ];\)\)\(\ \ \)\)\)

Text output - stability of rest points

\(\(\(
\)\(If[
        findrestpoints == 
          1, 

If[\ \ \(Dimensions[
                stblpnts]\)[\([1]\)] == 1, 
          Print["   Stable rest point:"], 
          Print["   Stable rest points:"]]\ ; 
If[\ \ \(Dimensions[stblpnts]\)[\([1]\)] == 0, Print["    None"], 
          Print[SequenceForm["   ", 
              ColumnForm[
                Table[PaddedForm[stblpnts[\([i]\)], {7, 7}, 
                    NumberPadding -> {" ", " "}], {i, 
                    1, \(Dimensions[
                        stblpnts]\)[\([1]\)]}]]]\ ]\ ]; 

If[\ \ \(Dimensions[unstbl]\)[\([1]\)] == 1, 
          Print["   Unstable rest point:"], 
          Print["   Unstable rest points:"]]\ ; 
If[\ \ \(Dimensions[unstbl]\)[\([1]\)] == 0, Print["    None"], 
          Print[SequenceForm["   ", 
              ColumnForm[
                Table[PaddedForm[unstbl[\([i]\)], {7, 7}, 
                    NumberPadding -> {" ", " "}], {i, 
                    1, \(Dimensions[
                        unstbl]\)[\([1]\)]}]]]\ ]]]\ ;\)\)\)

Text output - the contour function

\(\(\(\(\(If[
        includecontourplot == 1 || \ \ graphcontour3d == 1, 
        Print[StyleForm[
            ColumnForm[{"", "Graphics"}], "Section"]]]\ ;\)\(
\)
    \)

    \(If[includecontourplot == 1 || \ graphcontour3d == 1, 
        Print[SequenceForm[\ "   Name of Contour Function = ", \ 
            SimplexContourFunction]\ ]; \ 
If[
          showcontourformula == 1, 
          Print["   Definition of Contour Function = ", 
            SequenceForm["   ψ", \ {x, y} // 
                MatrixForm, \ " = ", \ 
              SimplexContourFunction[{x1, x2}, {y1, 
                  y2}]\ ]]]; 
Print[
          SequenceForm["   Minimum Value = ", \ \ N[
              MinAttained]\ ]]\ \ ; \ 
Print[
          SequenceForm["   Maximum Value = ", \ \ N[
              MaxAttained]\ ]]\ ]\ ;\)\)\(\ \)\)\)

Generating the vector field diagram

If[
        vectorfielddiagram == 
          1,  
Show[  If[
            drawrestpoints == 1, {convf, simplexcolorbarvf,  Tri, 
              vectorf1, vectorf2, DOTS, WDots}, {convf, simplexcolorbarvf,  
              Tri, vectorf1, vectorf2}], 
          DisplayFunction ->  $DisplayFunction, 
 
          AspectRatio ->  Automatic,  ImageSize ->  {500, 500}, 
          DefaultFont -> {"Palatino-Roman", 10}, 
          Epilog -> {Text[strategy[2], {- .025, 0}], 
              Text[strategy[1], {- .025, 1}], 
              Text[strategy[3], {0, 1.025}], 
              Text[strategy[4], {1, 1.025}]} ] ] ;

Generating the phase diagram on simplex

If[phasediagram == 1,  
        Show[{conph, Tri, e, ecycle, simplexcolorbarph, extralines, uddots, 
            udwdots, extratext, ar, arcycle, DOTS, WDots}, 
          DisplayFunction  ->  $DisplayFunction, 
          ImageSize ->  {500, 500}, 
          PlotRange -> {{- .08, 1.1}, {- .08, 1.1}}, 
          AspectRatio ->  Automatic, 
          DefaultFont -> {"Palatino-Roman", 10}, 
          Epilog -> {Text[strategy[2], {- .025, 0}], 
              Text[strategy[1], {- .025, 1}], 
              Text[strategy[3], {0, 1.025}], 
              Text[strategy[4], {1, 1.025}]} ] ] ;

Generating the 3D plot of the contour function

If[graphcontour3d ==  1,  
        Show[confunplt, 
          Which[findrestpoints == 1  &&  grayline == 1, {Tri3d, 
              SeveralLines, SeveralPoints, edgeline1, edgeline2, edgeline3, 
              edgeline4}, 
            findrestpoints == 1  &&  grayline == 0, {Tri3d, 
              SeveralLines, SeveralPoints}, 
            findrestpoints == 0  &&  grayline == 1, {Tri3d, 
              edgeline1, edgeline2, edgeline3, edgeline4}, 
            findrestpoints == 0  &&  grayline == 0, Tri3d], 
          Graphics3D[{Text[
                StyleForm[strategy[2], 
                  FontFamily -> "Palatino-Roman", 
                  FontSize -> 10], {- .04, - .04,  
                  MinAttained - distancefrommin}], 
Text[
                StyleForm[strategy[1], 
                  FontFamily -> "Palatino-Roman", 
                  FontSize -> 10], {- .04, 
                  1.04,   MinAttained - 
                    distancefrommin}], 
Text[
                StyleForm[strategy[3], 
                  FontFamily -> "Palatino-Roman", 
                  FontSize -> 10], { .04, 
                  1.06,   MinAttained - distancefrommin}], 
              Text[StyleForm[strategy[4], 
                  FontFamily -> "Palatino-Roman", 
                  FontSize -> 10], {1.04, 
                  1.06,   MinAttained - distancefrommin}]} ], 
          DisplayFunction  ->  $DisplayFunction,  
          ImageSize ->  {500, 500}]] ;

Print[
        SequenceForm["Runtime: ", 
          If[Quotient[TimeUsed[]  -  time1, 60] > 9, 
            Quotient[TimeUsed[]  -  time1, 60], 
            PaddedForm[Quotient[TimeUsed[]  -  time1, 60], 1, 
              NumberPadding ->  {"0", "0"}, 
              SignPadding ->  False]], ":", 
          If[Round[Mod[TimeUsed[]  -  time1, 60]] > 9, 
            Round[Mod[TimeUsed[]  -  time1, 60]], 
            PaddedForm[Round[Mod[TimeUsed[]  -  time1, 60]], 1, 
              NumberPadding ->  {"0", "0"}, 
              SignPadding ->  False]]]] ;

\(\(\(
\)\( (*\(\(--\(--\(--\(--\(--Returning\)\)\)\)\)\ to\ original\ output\ \(\(\(production--\)--\)--\)\)\(-\)*) \)\(
\)\(
\)\(\($DisplayFunction = 
        CellPrint[\ 
            Cell[GraphicsData["PostScript", \ 
                DisplayString[#]], "Graphics"]] &\ ;\)

    \(\($Post = 
        If[# =!= \ Null, 
            CellPrint[
              Cell[BoxData[
                  MakeBoxes[#, 
                    StandardForm]], "Output"]]] &\ ;\)\(
\)
    \)

    \(Unprotect[Print]\ ;\)

    \(Print[x_]\  := \ 
        CellPrint[
          Cell[BoxData[MakeBoxes[x, StandardForm]]]]\ ;\)

    \(Print[x__]\  := Map[Print, {x}]\ ;\)

    \(\(Protect[Print]\ ;\)\(
\)
    \)

    \(Unprotect[StylePrint]\ ;\)

    \(StylePrint[x_]\  := \ 
        CellPrint[
          Cell[BoxData[ToBoxes[x]], "Output"]]\ ;\)

    \(StylePrint[x_, rest___]\  := \ 
        CellPrint[Cell[BoxData[ToBoxes[x]], rest]]\ ;\)

    \(Protect[StylePrint];\)


    \)\( (*\(--\(--\(--\(--\(--\(--\(--\(--code\)\)\)\)\)\)\)\)\ \(\(\(\(\(ends--\)--\)--\)--\)--\) - \(--\(---\)\)*) \)\)\)

Initialization Cells; Procedures; Known Issues

Preliminary definitions and calculations

Remove["Global`*"] ;

<< "Graphics`Arrow`" ;

Calculation of restricted equilibria and Nash equilibria

Clear[NashEq] ;

    
NashEq::usage = "is a procedure that finds Nash Eq. of a game \
when payoff vector field, F, as an argument." ;

    
\(\(\(NashEq::notes = "finds Nash Equilibria by first finding rest points of Rep. Dyn. and eleminating the ones that are not Nash Eq.";\)\(
\)
    \)\)
\(\(NashEq[F1_, F2_] := 
        Module[\ {usex, yusey, i, j, nnR, TnnTR, nnsonuc, nnfin, RawRestRep, 
            RestRep, grd1, Nash0, Nash1}, 
          usex = Table[i\  .1 -  .1, \ {i, 11}]\ ; \ 
          yusey = Table[i\  .1 -  .1, \ {i, 11}]\ ; 
i = 
            1\ ; \ j = 
            1\ ; \ 
nax = {nax1, \ 
              1 - nax1}\ ; 
nay = {nay1, 
              1 - nay1}\ ; \ 

 (*\ 
            Finds\ rest\ of\ replicator\ by\ "grid search"\ where\ grids\ are\ given\ by\ usex\ and\ yusey\ *) 

nnR = 
            Flatten[Table[
                Table[FindRoot[{0 == \(Replicator[nax, nay, F1, 
                            F2]\)[[1]], 
                      0 == \(Replicator[nax, nay, F1, 
                            F2]\)[[3]]}, {nax1, usex[[i]]}, {nay1, 
                      yusey[[j]]}, \ 
                    MaxIterations\  -> 
                      1000], {j, \(Dimensions[
                        yusey]\)[[1]]}], {i, \(Dimensions[usex]\)[[1]]}], 
              1]; 

 (*\ 
            This\ function\ is\ defined\ to\ eleminate\ double\ countings\ in\ nnR\ caused\ by\ grid\ \(\(search\)\(.\)\)\ *) \ \ 

TnnTR[H_] := 
            DeleteCases[
              Table[If[\((Abs[\((nax1\  /. \ 
                                H[[1, 
                                  1]])\) - \ \((nax1\  /. \ H[[i, 1]])\)] >  .01\ \  || \ 
                      Abs[\((nay1\  /. \ \ H[[1, 
                                  2]])\) - \ \((nay1\  /. \ H[[i, 2]]\ )\)] >  .01\ )\)\ , 
                  H[[i]], 0], {i, 
                  2, \(Dimensions[
                      H]\)[[1]]}], 
              0]; 

 (*\ 
            Using\ TnnTR\ we\ eliminate\ double\ countings\ by\ following\ two\ lines\ *) \ \ 

nnsonuc = 
            DeleteCases[
              FixedPointList[TnnTR, nnR], {}]\ ; 
nnfin = 
            Chop[Table[\(nnsonuc[[i]]\)[[1]], {i, \(Dimensions[
                      nnsonuc]\)[[1]]}], 
              10\^\(-6\)]; 

 (*\ 
            Following\ two\ lines\ eleminates\ the\ values\ outside\ the\ \(\(simplex\)\(.\)\)\ *) 

RawRestRep = 
            Chop[{{nax1, 1 - nax1}, {nay1, 1 - nay1}}\  /. \ nnfin, 
              10\^\(-6\)]; 
RestRep = 
            DeleteCases[
              Table[If[
                  Min[RawRestRep[[i]]] < \(-10\^\(-6\)\), outside, 
                  RawRestRep[[i]]], {i, \(Dimensions[RawRestRep]\)[[1]]}], 
              outside]; 

grd1 = {{1, 
                0}, {0, 1}}; 

 (*\ 
            Following\ three\ lines\ eleminate\ the\ rest\ point\ that\ are\ not\ Nash\ \(\(Eq\)\(.\)\)\ *) 

Nash0 = 
            Table[Table[
                If[\((F1[RestRep[\([j, 1]\)], RestRep[\([j, 2]\)]] . 
                            RestRep[\([j, 
                                1]\)] +  .00001\  ≥ \ \ F1[
                            RestRep[\([j, 1]\)], RestRep[\([j, 2]\)]] . 
                          grd1[\([i]\)]\  && \ 
                      F2[RestRep[\([j, 1]\)], RestRep[\([j, 2]\)]] . 
                            RestRep[\([j, 
                                2]\)] +  .00001\  ≥ \ \ F2[
                            RestRep[\([j, 1]\)], RestRep[\([j, 2]\)]] . 
                          grd1[\([i]\)])\), \ RestRep[\([j]\)], 
                  sil], {i, \(Dimensions[grd1]\)[\([1]\)]}], {j, \(Dimensions[
                    RestRep]\)[\([1]\)]}]\ ; 
Nash1 = 
            DeleteCases[
              Table[If[MatrixQ[Flatten[Nash0[\([i]\)], 1]], Nash0[\([i]\)], 
                  sil], {i, \(Dimensions[Nash0]\)[\([1]\)]}], 
              sil]\ ; 


Nash2 = 
            Table[Nash1[\([i, 
                  1]\)], {i, \(Dimensions[
                    Nash1]\)[\([1]\)]}]; \ \ Nash2\ ]\ ;\)\)

Calculation of rest points (smooth dynamics)

Clear[RestPointsSmooth]  ;

    
RestPointsSmooth::usage = "Calculates restpoint of smooth \
dynamics "  ;

    
RestPointsSmooth[Dynamic_, F1_, F2_,  RestSearchInit_,  
          TolLevel_] :=  
        Module[{sonuc, vfield, TR, fin, RAWresult},  
          srex = {srex1, 1 - srex1} ;  
          srey = {srey1, 1 - srey1} ;   
vfield[Dynamic, 
              x_, y_, k1_, k2_]  :=  
            Dynamic[x, y, k1, 
              k2] ; 

 (* 
            Finds restpoints by grid search *)  
          R = Table[
              FindRoot[{0 == vfield[Dynamic, srex, srey, F1, 
                        F2][[1]], 
                  0 == vfield[Dynamic, srex, srey, F1, 
                        F2][[3]]}, {srex1, RestSearchInit[[i, 1]]}, {srey1, 
                  RestSearchInit[[i, 
                    2]]},  
                MaxIterations  -> 
                  1000], {i, Dimensions[
                    RestSearchInit][[1]]}] ; 

 (* 
            Eleminates double countings *) 
TR[H_] := 
            DeleteCases[
              Table[If[(Abs[(srex1  /.  
                                H[[1, 
                                  1]]) -  (srex1  /.  H[[i, 1]])] >  .01   ||  
                      Abs[(srey1  /.   H[[1, 
                                  2]]) -  (srey1  /.  H[[i, 2]] )] >  .01 ) , 
                  H[[i]], 0], {i, 
                  2, Dimensions[
                      H][[1]]}], 
              0]; 
sonuc = 
            DeleteCases[FixedPointList[TR, R], {}]; 
fin = 
            Chop[Table[sonuc[[
                    i]][[1]], {i, Dimensions[
                      sonuc][[1]]}], 
              (10^(-6))] ; 

RAWresult = Chop[{{srex1, 1 - srex1}, {srey1, 1 - srey1}}  /.  fin, (10^(-6))] ;  
          midresult = 
            DeleteCases[
              Table[If[
                  Min[RAWresult[[i]]] < -(10^(-6)), outside, 
                  RAWresult[[i]]], {i, Dimensions[RAWresult][[1]]}], 
              outside] ;  
          result = 
            DeleteCases[
              Table[If[
                  Max[Abs[
                        Dynamic[
                          midresult[[i, 
                            1]], 
                          midresult[[i, 
                            2]], F1, F2]]]  >  
                    TolLevel,  notrest, 
                  midresult[[i]]], {i, Dimensions[midresult][[1]]}], 
              notrest] ; 
  result ] ;

Specification of rest points

Clear[RestPoints] ;
\(\(\(RestPoints::usage = "Calculation of rest points of dynamics"\ ;\)\(
\)
    \)\)
RestPoints[Dynamic_, F1_, F2_, RestSearchInit_, TolLevel_] := 
        Which[RPCharacterization[Dynamic] ===  Automatic, 
          RestPointsSmooth[Dynamic, F1, F2,  RestSearchInit, TolLevel],  
          RPCharacterization[Dynamic] ===  Nash,  NashEq[F1, F2],  
          RPCharacterization[Dynamic] ===  RestrictedEquilibria,  
          RestPointsSmooth[Replicator, F1, F2,  RestSearchInit,  
            TolLevel]] ;

Testing stability of rest points

Clear[
        TestStability] ;
\(TestStability::usage = "tests stability of rest points found by Respoints[Dynamic,F] function."\ ;\)
TestStability::notes = "Its arguments are name of the \
dynamic, payoff vector field, final time which will be used to see whether \
solution trajectory is outside neighborhood of the rest point. This \
neighborhood is defined by delta. eps defines another neighborhood of rest \
point from which an initial condition drawn. numtest gives the number of \
initial points for which tests will be conducted." ;

TestStability[Dynamic_, F1_, F2_, bigT_,  numtest_,  
          eps_,   delta_, RestSearchInit_, TolLevel_]:=Module[{xtest0, vfield, G1, G3, dim, Stability}, testrest = 
                RestPoints[Dynamic, F1, F2, RestSearchInit, TolLevel] ;   
Table[
                For[j = 1, 
                  j < numtest + 1, j++, (Label[AA]; θ =  
                      Random[Real, {0, 2 Pi}]; 
                    xtest0[i, 
                        j] = {testrest[[i, 1, 
                            1]], 
                          testrest[[i, 2, 
                            1]]} +  
                        eps {Cos[θ], Sin[θ]};  
                    If[Min[xtest0[i, j]] < 0 || Max[xtest0[i, j]] > 1,  
                      Goto[AA]])], {i, Dimensions[
                      testrest][[1]]}];

xtest = 
                Table[xtest0[i, 
                    j], {i, Dimensions[
                        testrest][[1]]}, {j, numtest}];

Clear[vfield];
vfield[Dynamic,  x_, y_, m1_, 
                  m2_]  :=  Dynamic[x, y, m1, m2] ;

  
trex = {trex1[t], 1 - trex1[t]}; trey = {trey1[t], 1 - trey1[t]};

G1=Table[NDSolve[{trex1'[t]==vfield[Dynamic, trex, trey, F1, 
                                  F2][[1]],trey1'[t]==vfield[Dynamic, trex, trey, F1, 
                                  F2][[3]],trex1[
                                0] == xtest[[i, 
                                  j]][[1]],trey1[
                                0] == xtest[[i, 
                                  j]][[2]]},{trex1, trey1},{t, 0, bigT}],{i, Dimensions[
                          testrest][[1]]},{j, Dimensions[
                          xtest][[2]]}] ;


  

G3 = {{trex1[bigT], 
                      1 - trex1[bigT]}, {trey1[bigT], 
                      1 - trey1[bigT]}}  /.  G1 ;

Dimensions[G3];
Stability=Table[If[((∑_(i = 1)^2) (Flatten[
                                        testrest[[j]]]-Flatten[G3[[j, k, 
                                        1]]])[[2 i - 1]]^2)^(1/2) < delta,1,0],{j, Dimensions[
                          G3][[1]]},{k, Dimensions[
                          G3][[2]]}];

dim = Dimensions[
                    Stability][[2]];

  
stblpnts=DeleteCases[Table[If[ (∑_(j = 1)^dim) Stability[[i, 
                                    j]]==dim,testrest[[i]] ],{i, Dimensions[Stability][[1]]}],Null] ;
unstbl=DeleteCases[Table[If[ (∑_(j = 1)^dim) Stability[[i, 
                                    j]]≠dim,testrest[[i]] ],{i, Dimensions[Stability][[1]]}],Null] ; 
{stblpnts, unstbl}] ;

Known issues

Rest Points of BR dynamics and Logit Dynamics

Strictly speaking  there is no built-in best response dynamic in this program, but we can approximate it using Logit[η] dynamic with a small value of  η.  But when η is very small, the FindRoot command (which is used to find the rest points, via Newton's method) does not work properly, and hence the rest points of approximate BR dynamic cannot be found accurately.  If you are using this dynamic, and "tuning" the rest point search does not help you, you may want to add rest point manually using the "Extra graphic elements" portion of the "Specification of output - Phase diagram" section.

Lyapunov functions for Stable Games

The Lyapunov functions for the replicator and logit dynamics in stable games may still need some work.

Credits

Dynamo 2x2
Version 0.1.5  (April 11, 2007)
by William H. Sandholm and Emin Dokumaci

mail:			Department of Economics
			University of Wisconsin
			1180 Observatory Drive
			Madison WI 53706

e-mail:			whs@ssc.wisc.edu
			edokumaci@wisc.edu
	
websites:		http://www.ssc.wisc.edu/~whs/
			http://www.ssc.wisc.edu/~edokumac/
		
Dynamo website:		http://www.ssc.wisc.edu/~whs/dynamo
    
Financial support from the National Science Foundation under Grants SES-0092145 and SES-0617753 is gratefully acknowledged.