/* Maxima is a GNU free software package for algebraic computations.	  */
/* For downloads and information go to:	 	*/
/*        http://maxima-project.org/		*/
/* The following Maxima code was developed by Tyler Smith. 	*/

/* Reduce the precision of floating point values when they are printed (full precision is used in calculations) */
fpprintprec : 3;

/* Section 7.2: vector fields */
load(plotdf);
eq1: r1 * n1 ;
eq2: r2 * n2 ;

/* Figure 7.1 */
plotdf([eq1, eq2], [n1, n2], [n1, -1, 1], [n2, -1, 1], [parameters, "r1=1,r2=2"]);
plotdf([eq1, eq2], [n1, n2], [n1, -1, 1], [n2, -1, 1], [parameters, "r1=1,r2=-2"]);
plotdf([eq1, eq2], [n1, n2], [n1, -1, 1], [n2, -1, 1], [parameters, "r1=-1,r2=-2"]);

/* The plots can also be manipulated with the _config_ dialog, followed by _replot_ in the plot window */

/* Section 7.3: equilibria and null clines */
kill(all);

/* Figure 7.2 */
eq1: a * n1 + b * n2 ;
eq2: c * n1 + d * n2 ;

/* System equilibrium */
/* Both solve, algsys and linsolve work, I don't know if one is better than the others... */
algsys([eq1, eq2], [n1, n2]);
linsolve([eq1, eq2], [n1, n2]);
solve([eq1, eq2], [n1, n2]);

/* Null clines */
nc1: linsolve([eq1], [n2]);
nc2: linsolve([eq2], [n2]);

a : 1; b : 2; c : -1; d : 3;
plot2d([rhs(nc1[1]), rhs(nc2[1])], [n1, -1, 1], [y, -1, 1], [ylabel, n2], [legend, "n1 null cline", "n2 null cline"]);

/* Section 7.3.2: eigenvalues, eigenvectors, and vector fields */
/* Figure 7.3a */

a : 2; b : 1; c : 1/2; d : 1;

load("eigen");
M : coefmatrix([eq1, eq2], [n1, n2]);
eigM : eivects(M);

/* The eigenvalues are found in the first sublist returned by eivects. Extraction is a little convoluted... first()
extracts list elements. ev() evalutes the expression using the numerical values of all the parameters. float() converts
the result to a decimal number */
ev1: float(ev(first(first(first(eigM)))));
ev2: float(ev(second(first(first(eigM)))));

/* eivects returns the eigenvectors with the x value scaled to one, so the slope is the second element of the list */
plot2d([first(rest(eigM))[2] * x, second(rest(eigM))[2] * x], [x, -1, 1], [y, -1, 1], [legend, concat("r = ", ev1),\
  concat("r = ", ev2)]);

/* A bug(?) in plotdf prevents me from plotting the eigenvectors directly on the vector field, so we have to use a
separate plot for the eigenvectors and the vector field */

plotdf([eq1, eq2], [n1, n2], [n1, -1, 1], [n2, -1, 1]);

/* Update the parameters and replot the figures */
/* Figure 7.3b */

a : -2; b : -1; c : 1/2; d : 1;
ev1: float(ev(first(first(first(eigM)))));
ev2: float(ev(second(first(first(eigM)))));

plot2d([first(rest(eigM))[2] * x, second(rest(eigM))[2] * x], [x, -1, 1], [y, -1, 1], [legend, concat("r = ", ev1),\
  concat("r = ", ev2)]);

plotdf([eq1, eq2], [n1, n2], [n1, -1, 1], [n2, -1, 1]);

/* Figure 7.3c */

a : -2; b : -1; c : -1/2; d : -1;
ev1: float(ev(first(first(first(eigM)))));
ev2: float(ev(second(first(first(eigM)))));

plot2d([first(rest(eigM))[2] * x, second(rest(eigM))[2] * x], [x, -1, 1], [y, -1, 1], [legend, concat("r = ", ev1),\
  concat("r = ", ev2)]);

plotdf([eq1, eq2], [n1, n2], [n1, -1, 1], [n2, -1, 1]);

/* Example: Metastasis of Malignant Tumors */

dC : -d1 * C - b * C;
dI : b * C - d2 * I + p * I;

M : coefmatrix([dC, dI], [C, I]);
Mevals : eivals(M)[1];

/* Figure 7.6 */
eq1: a * n1 + b * n2 ;
eq2: c * n1 + d * n2 ;

a : -2; b : -1; c : 1; d : -2;
plotdf([eq1, eq2], [n1, n2], [n1, -1, 1], [n2, -1, 1]);

a : 2; b : -1; c : 1; d : 2;
plotdf([eq1, eq2], [n1, n2], [n1, -1, 1], [n2, -1, 1]);

/* Example: evolution by sexual selection */
/* 7.17 */
dz : Gz * a * p - Gz * c * z;
dp : -Gp * b * p;
M : coefmatrix([dz, dp], [z, p]);
Mvals : eivals(M)[1];

/* 7.19 */
dz : G_z *( a * p - c * z) - B * b * p;
dp : B * (a * p - c * z) - G_p * b * p;
M : coefmatrix([dz, dp], [z, p]);
Mvals : eivals(M)[1];

a : 0.8; b : 0.5; c : 0.9; B : 0.2; G_z : 0.5; G_p : 0.4;
plotdf([dz, dp], [z, p], [z, -10, 10], [p, -10, 10], [trajectory_at, -7.55, 7.55], [nsteps, 500]);
/* Plot vs time using the GUI interface */

a : 5; b : 0.5; c : 0.9; B : 0.2; G_z : 0.5; G_p : 0.4;
plotdf([dz, dp], [z, p], [z, -10, 10], [p, -10, 10], [trajectory_at, 0.1, -0.1], [nsteps, 500]);

/* To quit: */

quit();