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

User-defined Parameters

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

Off[General::spell1];

Choice of game

Choose a normal form game here.

M1=(({{{0, 0, 0}, {1, 1, 0}}, {{1, 0, 1}, {0, 1, 1}}})) ; M2=(({{{0, 1, 1}, {1, 0, 1}}, {{1, 1, 0}, {0, 0, 0}}})) ;

Nonlinear population games

To choose a nonlinear game, set nonlineargame = 1 and define F1[x_,y_,z_] , F2[x_,y_,z_] and F3[x_,y_,z_] as you like. In doing so, the components of x should be entered as x[[1]], x[[2]],  and the components of y and z should be entered similarly.  The definitions of A[. , ] are used in the normal form case.

nonlineargame  = 0 ;

A[1, 1]=(({{M1[[1, 1, 
                  1]], M1[[1, 2, 
                  1]]}, {M2[[1, 1, 
                  1]], M2[[1, 2, 
                  1]]}})) ;
A[1, 2]=(({{M1[[2, 1, 
                  1]], M1[[2, 2, 
                  1]]}, {M2[[2, 1, 
                  1]], M2[[2, 2, 
                  1]]}})) ;
A[2, 1]=(({{M1[[1, 1, 
                  2]], M1[[2, 1, 
                  2]]}, {M2[[1, 1, 
                  2]], M2[[2, 1, 
                  2]]}}));
A[2, 2]=(({{M1[[1, 2, 
                  2]], M1[[2, 2, 
                  2]]}, {M2[[1, 2, 
                  2]], M2[[2, 2, 
                  2]]}}));
A[3, 1]=(({{M1[[1, 1, 
                  3]], M1[[1, 2, 
                  3]]}, {M1[[2, 1, 
                  3]], M1[[2, 2, 
                  3]]}}));
A[3, 2]=(({{M2[[1, 1, 
                    3]], M2[[1, 2, 
                    3]]}, {M2[[2, 1, 
                    3]], M2[[2, 2, 
                    3]]}}));

F1[x_, y_, 
        z_]  :=  {z . (A[1, 1] . y ), z . (A[1, 2] . y )};  
    F2[x_, y_, z_]  := {z . (A[2, 1] . x ), 
        z . (A[2, 2] . x )} ;
F3[x_, y_, 
        z_] :=  {x . (A[3, 1] . y ), x . (A[3, 2] . y )} ;  
    F[x_, y_, z_] := {F1[x, y, z], F2[x, y, z], F3[x, y, z]} ;

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" ;
strategy[5] = "I" ;
strategy[6] = "O" ;

Payoff-related definitions

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

The population's average payoff

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

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

    
Fbar[3][x_, y_, z_, 
          F3_] :=  (∑_(i = 1)^2) z[[
              i]] F3[x, y, z][[
              i]] ;

Excess payoff vector

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

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

    
Fhat[3][x_, y_, z_, F3_] := 
        F3[x, y, z]  -  Fbar[3][x, y, z, 
            F3] ;

Vector of positive parts of excess payoffs

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

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

    
Fhatplus[3][x_, y_, z_, 
          F3_] :=  {Max[
            0, Fhat[3][x, y, z, 
                F3][[1]]], 
          Max[0, Fhat[3][x, y, z, 
                F3][[2]]]} ;

Projected payoff vector

PhiF[x_, y_, z_, F1_, F2_, F3_] := 
        Flatten[{F1[x, y, 
                z] - {1/2, 
                  1/2} (∑_(i = 1)^2) F1[x, y, z][[i]], 
            F2[x, y, 
                z] - {1/2, 
                  1/2} (∑_(i = 1)^2) F2[x, y, z][[i]], 
            F3[x, y, 
                z] - {1/2, 
                  1/2} (∑_(i = 1)^2) F3[x, y, z][[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. The definition of the projection dynamic is only correct on the interior of the simplex, so boundary trajectories of the projected payoff dynamic must be added by hand.   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 ;

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.  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;

    
Combined[x_, y_, z_, F1_, F2_, F3_] := 
        weight[1] dyn[1][x, y, z, F1, F2, F3] + 
          weight[2] dyn[2][x, y, z, F1, F2, F3];

Definitions of dynamics

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

Replicator Dynamic

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

RPCharacterization[Replicator] = 
        Automatic ;

Maynard Smith Replicator Dynamic

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

RPCharacterization[MSReplicator] = 
        Automatic ;

Logit Dynamic

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

RPCharacterization[Logit[_]] = 
        Automatic ;

I-Logit Dynamic

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

RPCharacterization[ILogit[_]] = 
        Automatic ;

BNN dynamic

BNN[x_, y_, z_, F1_, F2_, 
          F3_] := {Fhatplus[1][x, y, z, 
                F1][[1]] - 
            x[[1]] (∑_(i = 1)^2) Fhatplus[1][x, y, z, F1][[
                  i]], Fhatplus[1][x, y, z, 
                F1][[2]] - 
            x[[2]] (∑_(i = 1)^2) Fhatplus[1][x, y, z, F1][[
                  i]], Fhatplus[2][x, y, z, 
                F2][[1]] - 
            y[[1]] (∑_(i = 1)^2) Fhatplus[2][x, y, z, F2][[
                  i]], Fhatplus[2][x, y, z, 
                F2][[2]] - 
            y[[2]] (∑_(i = 1)^2) Fhatplus[2][x, y, z, F2][[
                  i]], Fhatplus[3][x, y, z, 
                F3][[1]] - 
            z[[1]] (∑_(i = 1)^2) Fhatplus[3][x, y, z, F3][[
                  i]], Fhatplus[3][x, y, z, 
                F3][[2]] - 
            z[[2]] (∑_(i = 1)^2) Fhatplus[3][x, y, z, F3][[
                  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,z,Fi], which is defined as Fhat[i][x,y,z,Fi] = Fi[x,y,z] - Fbar[i][x,y,z,Fi], for each population i = 1, 2,3.  (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_, z_, F1_, F2_, 
          F3_] :=  {(Max[
                0, Fhat[1][x, y, z, 
                    F1][[1]]])^2, 
(Max[
                0, Fhat[1][x, y, z, 
                    F1][[2]]])^2, (Max[0, Fhat[2][x, y, z, 
                    F2][[1]]])^2, (Max[0, Fhat[2][x, y, z, 
                    F2][[2]]])^2, (Max[0, Fhat[3][x, y, z, 
                    F3][[1]]])^2, (Max[0, Fhat[3][x, y, z, 
                    F3][[2]]])^2};

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

RPCharacterization[EP] = Nash  ;

Pairwise Difference Dynamic

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

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

RPCharacterization[PD] = Nash ;

Pairwise Comparison Dynamic

Here one specifies the vector field rho1 that defines an excess payoff dynamic.

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

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

RPCharacterization[PC] = Nash ;

Projection Dynamic

closetozero = 10^(-7) ;

Projection[x_, y_, z_, F1_, F2_, 
          F3_] := {Evaluate[
            If[x[[1]] > 
                  closetozero  ||  (F1[x, y, 
                        z][[1]] ≥ F1[x, y, z][[2]]), 
              If[x[[2]] > 
                    closetozero  ||  (F1[x, y, 
                          z][[2]] ≥  F1[x, y, 
                          z][[1]]), F1[x, y, 
                      z][[1]] - ((1/2) ((∑_(j = 1)^2) F1[x, y, z][[
                            j]])), 0], 0]], 
          Evaluate[
            If[x[[2]] > 
                  closetozero  ||  (F1[x, y, 
                        z][[2]] ≥  F1[x, y, 
                        z][[1]]), 
              If[x[[1]] > 
                    closetozero  ||  (F1[x, y, 
                          z][[1]] ≥  F1[x, y, 
                          z][[2]]), F1[x, y, 
                      z][[2]] - ((1/2) ((∑_(j = 1)^2) F1[x, y, z][[
                            j]])), 0], 0]], 
          Evaluate[
            If[(y[[1]] > 
                    closetozero  ||  (F2[x, y, 
                          z][[1]] ≥ F2[x, y, 
                          z][[2]])) , If[(y[[2]] > 
                      closetozero  ||  (F2[x, y, 
                            z][[2]] ≥  F2[x, y, 
                            z][[1]])), F2[x, y, 
                      z][[1]] - ((1/2) ((∑_(j = 1)^2) F2[x, y, z][[
                            j]])), 0], 0]], 
          Evaluate[
            If[(y[[2]] > 
                    closetozero  ||  (F2[x, y, 
                          z][[2]] ≥  F2[x, y, 
                          z][[1]])), If[(y[[1]] > 
                      closetozero  ||  (F2[x, y, 
                            z][[1]] ≥ F2[x, y, 
                            z][[2]])) , F2[x, y, 
                      z][[2]] - ((1/2) ((∑_(j = 1)^2) F2[x, y, 
                              z][[j]])), 0], 0]], 
          Evaluate[
            If[(z[[1]] > 
                    closetozero  ||  (F3[x, y, 
                          z][[1]] ≥ F3[x, y, 
                          z][[2]])) , If[(z[[2]] > 
                      closetozero  ||  (F3[x, y, 
                            z][[2]] ≥  F3[x, y, 
                            z][[1]])), F3[x, y, 
                      z][[1]] - ((1/2) ((∑_(j = 1)^2) F3[x, y, 
                              z][[j]])), 0], 0]], 
          Evaluate[
            If[(z[[2]] > 
                    closetozero  ||  (F3[x, y, 
                          z][[2]] ≥  F3[x, y, 
                          z][[1]])), If[(z[[1]] > 
                      closetozero  ||  (F3[x, y, 
                            z][[1]] ≥ F3[x, y, 
                            z][[2]])) , F3[x, y, 
                      z][[2]] - ((1/2) ((∑_(j = 1)^2) F3[x, y, 
                              z][[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_, z_, F1_, F2_, 
            F3_] := {} ;              ;

RPCharacterization[Other] = Automatic  ;

Specification of solution trajectories

Instructions for specification of solution trajectories

data specfies the solution trajectories to be drawn. ex,: data= {{{ .989,  .99,  .99}, 20, {Thickness[.004], Hue[1]}},
  {{ .9905,  .99,  .99}, 20, {Thickness[.004] ,Hue[.65]}}}. The first entry in the example, {{ .989,  .99,  .99}, 20, {Thickness[.004], Hue[1]}} is for the first solution trajectory.  The entry { .989,  .99,  .99} describes the initial condition of the solution.  Each initial condition is entered as a list {x2, y2, z2} containing the weight placed on each population's second strategy.   The entry 20 determines the time span over which which the differential equational will be solved, [0, 20].  Hue[1] determines in which color the solution trajectory will be drawn, in this case red. Any color command can be used, i.e. GrayLevel, Hue, RGBColor. This structure must be followed in order for the program run correctly (however you can take out or add graphics primitives like Dashing[{...}], Thickness  in the last entry of each sol trajectory. i.e Hue[1], or {Thickness[.004], Hue[1]}, {Dashing[{.01,.03}], Thickness[.004], Hue[1]}, etc). If only one solution trajectory is needed, do not forget the outermost brackets, i.e. {  {{.011, .989, .01, .99, .01, .99}, 20, {Thickness[.004], Hue[1]}} }. 
Dotdata specifies dots to be drawn at rest points, etc. ex.: dotdata={{.01, {.7, .2, .3}, GrayLevel[.2]}, {...........}}. creates a spehere with radius .0, centered at {.3, .8, .7} and with color GrayLevel[.2 ].  The color that you will see is influenced by a number of factors.  (Simulated illumination determines the shading of polygons in three-dimensional pictures. The shading of a particular polygon is computed as a sum of contributions from point light sources, plus a contribution from ambient light. ) 
ppoints determines the number of line segments used when drawing the solution trajectories.  Increase its value if the trajectories have kinks.

data =  {{{ .5,  .51,  .5}, 
            24, {Thickness[ .005], Hue[1]}}, {{ .5,  .5,  .51}, 
            9, {Thickness[ .005], Hue[ .7]}}} ;

    
dotdata = {{ .01, { .5,  .5,  .5}, 
            GrayLevel[ .2]} };

    
ppoints = 1000 ;

Specify the value of frontedge to determine whether edge LI or edge RI is in the foreground of the output.  This choice determines the viewpoint used by Mathematica (LIviewpoint or RIviewpoint) as well as which edges of the cube are dashed.

({{LI , RI }})

frontedge = RI  ;

    
LIviewpoint = {1.558, 2.877, 0.862} ;
RIviewpoint = {-1.558, 2.877, 0.862} ;

Program Core

Solving the differential equation

time1 = TimeUsed[] ;

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

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

Plotting the solutions

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

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

Drawing Triangle and Dots

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

<< Graphics`Shapes`

simplex = 
        Which[frontedge ===  LI, 
          Graphics3D[{Line[{{0, 1, 0}, {1, 1, 0}}], 
              Line[{{1, 1, 0}, {1, 1, 1}}], Line[{{1, 1, 1}, {0, 1, 1}}], 
              Line[{{0, 1, 1}, {0, 1, 0}}], Line[{{1, 1, 1}, {1, 0, 1}}], 
              Line[{{1, 0, 1}, {0, 0, 1}}], Line[{{0, 0, 1}, {0, 1, 1}}], 
              Line[{{1, 1, 0}, {1, 0, 0}}], 
              Line[{{1, 0, 0}, {1, 0, 1}}], {Dashing[{ .02,  .01}], 
                Line[{{0, 0, 0}, {1, 0, 0}}], Line[{{0, 0, 0}, {0, 1, 0}}], 
                Line[{{0, 0, 0}, {0, 0, 1}}]}}], frontedge ===  RI,  
          Graphics3D[{Line[{{0, 1, 0}, {1, 1, 0}}], 
              Line[{{1, 1, 0}, {1, 1, 1}}], Line[{{1, 1, 1}, {0, 1, 1}}], 
              Line[{{0, 1, 1}, {0, 1, 0}}], Line[{{1, 1, 1}, {1, 0, 1}}], 
              Line[{{1, 0, 1}, {0, 0, 1}}], 
              Line[{{0, 0, 1}, {0, 1, 1}}], {Dashing[{ .02,  .01}], 
                Line[{{1, 1, 0}, {1, 0, 0}}], Line[{{1, 0, 0}, {1, 0, 1}}], 
                Line[{{0, 0, 0}, {1, 0, 0}}]}, Line[{{0, 0, 0}, {0, 1, 0}}], 
              Line[{{0, 0, 0}, {0, 0, 1}}]}]] ;

    
kureler = 
        Table[Graphics3D[{SurfaceColor[ dotdata[[i, 3]] ], EdgeForm[], 
              TranslateShape[
                Sphere[dotdata[[i, 1]], 20, 
                  15],  {1 - 
                    dotdata[[i, 2, 2]], 
                  1 - dotdata[[i, 2, 
                      3]], 
                  1 - dotdata[[i, 2, 
                      1]]} ] }] , {i, Dimensions[
                dotdata][[1]]}] ;

Drawing Solutions to Diff Eqs if there is no contour plot

e = 
        Table[ParametricPlot3D[
            Evaluate[
              Join[{y1[t], z1[t], 
                    x1[t]} /. G[[i, 
                    1]] , { 
                  data[[i, 3]] } ]  ], {t, 0, data[[i, 2]]},  
            DisplayFunction  ->  Identity, 
            PlotPoints ->  ppoints, Boxed ->  False,  
            Axes -> False, 
            PlotRange -> {{- .08, 1.08}, {- .08, 
                  1.08}, {- .08, 1.08}}, 
            AspectRatio -> 
              Automatic], {i, Dimensions[
                data][[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[
                          5], ""}, {"", ""}, {strategy[3], 
                      strategy[4]}}, M1]] // DisplayForm, 
            GridBox[{{"    "}, {"    "}, {"    "}, {"    "}}]\  // DisplayForm\ , 
            GridBox[{{""}, {""}, {""}, {strategy[1]}, {strategy[
                      2]}}]\  // DisplayForm\ , 
            GridBox[Join[{{"            "\ strategy[
                          6], ""}, {"", ""}, {strategy[3], 
                      strategy[4]}}, M2]] // 
              DisplayForm]]]\ ;\)\(
\)
    \)\)
Print[
        SequenceForm["   Payoff vector field:  ", "F", {x, y, z} // 
            MatrixForm,  " = ", 
          Flatten[F[{x1, x2}, {y1, y2}, {z1, z2}]]  //  
            MatrixForm]] ;

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, z} // 
            MatrixForm, \ " = ", 
          Simplify[Dynamic[{x1, x2}, {y1, y2}, {z1, z2}, F1, F2, F3]] // 
            MatrixForm\ ]\ ];\)\)\(\ \ \)\)\)

Generating the phase diagram on simplex

Show[e, kureler, simplex, 
        Which[frontedge ===  LI, 
          Graphics3D[{Text[strategy[1], {- .02, 1.02,  .98}], 
              Text[strategy[2], {- .03, 1.02,  .02}], 
              Text[strategy[4], { .01, 1.02, - .04}], 
              Text[strategy[3], {1, 1.02, - .04}], 
              Text[strategy[5], {0,  .98, 1.04}], 
              Text[strategy[6], {0,  .08, 1.02}]}], frontedge ===  RI,  
          Graphics3D[{Text[strategy[1], {1.02, 1.03,  .98}], 
              Text[strategy[2], {1.02, 1.02,  .02}], 
              Text[strategy[4], { .01, 1.02, - .04}], 
              Text[strategy[3], { .98, 1.02, - .04}], 
              Text[strategy[5], {1,  .98, 1.04}], 
              Text[strategy[6], {1, 0.08, 1.03}]}]], 
        PlotRange ->  {{- .08, 1.08}, {- .08, 
              1.08}, {- .08, 1.08}}, 
        DisplayFunction  ->  $DisplayFunction, 
        AspectRatio ->  Automatic, 
        DefaultFont -> {"Palatino-Roman", 10},  
        ViewPoint ->  
          Which[frontedge ===  LI,  LIviewpoint,  frontedge ===  RI,  
            RIviewpoint], AmbientLight ->  Hue[ .55], 
        LightSources ->  {{{ .5, 1, 1}, 
              RGBColor[1, 0, 0]}, {{1,  .8,  .5}, 
              RGBColor[0, 1, 0]}, {{1, 1, 1}, 
              RGBColor[0, 0, 1]}}] ;

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--\)--\)--\)--\)--\) - \(--\(---\)\)*) \)\)\)

Credits

Dynamo 2x2x2
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.