# 4-dimensional Lotka-Volterra plot

I was referring to the question I posted earlier in Mathematics forum, but then I came to see this example on Wikipedia. How did they plot the system in phase space?

The closet relevant topic I have done is the Lorenz Equation, where I used NDSolve and then ParametricPlot3D evaluated at the fix points.

• It would be very helpful if you can at least include the setup of the equations you've setup. – Pillsy Apr 17 '17 at 19:16
• @Pillsy Hi, I am trying to model $\frac{dx_i}{dt}=r_ix_i\left(1-\sum_{j=1}^{N}a_{ij}x_j\right)$, but I have only been using Mathematica more than a month. What commend should I be using? – Vanya Apr 17 '17 at 19:39
• What initial conditions are you using? – Pillsy Apr 17 '17 at 19:40
• @Pillsy Not clear yet. I will be messing around with the parameters as well as the initial conditions to try to simulate the same result in here. But I think I may know the system more if I start from lower dimensions. – Vanya Apr 17 '17 at 19:43
• For future questions, try to include some code so that people don't have to type all the equations, constants, etc. to provide help. Also, folks like to see that the poster has made some efforts themselves. – Chris K Apr 18 '17 at 1:52

Here's how to numerically solve the model.

nsp = 4;

{r[1], r[2], r[3], r[4]} = {1, 0.72, 1.53, 1.27};
{k[1], k[2], k[3], k[4]} = {1, 1, 1, 1};
amat = {{1, 1.09, 1.52, 0}, {0, 1, 0.44, 1.36}, {2.33, 0, 1, 0.47}, {1.21, 0.51, 0.35, 1}};
Do[a[i, j] = amat[[i, j]], {i, nsp}, {j, nsp}];

eqns = Table[x[i]'[t] ==
r[i]*x[i][t]*(1 - Sum[a[i, j]*x[j][t]/k[i], {j, nsp}]), {i, nsp}];
ics = Table[x[i][0] == 0.1, {i, nsp}];
unks = Table[x[i], {i, nsp}];

tmax = 10000;
sol = NDSolve[{eqns, ics}, unks, {t, 0, tmax}][[1]];

ParametricPlot3D[Evaluate[{x[1][t], x[2][t], x[3][t]} /. sol], {t, 100, tmax},
AxesLabel -> {"x1", "x2", "x3"}, PlotPoints -> 200]


I don't know how to color the ParametricPlot3D according to x[4][t] -- that could be an interesting question if it hasn't already been answered on the site.

Update

Here's how to color according to x[4][t], using the slot #4 which represents t here:

ParametricPlot3D[Evaluate[{x[1][t], x[2][t], x[3][t]} /. sol], {t, 100, tmax},
AxesLabel -> {"x1", "x2", "x3"}, PlotPoints -> 400,
ColorFunction -> (ColorData["SunsetColors", x[4][#4] /. sol] &),
ColorFunctionScaling -> False]


• I'm not particularly pleased with how I set up the a[i,j] terms here. Should be a more elegant way... – Chris K Apr 18 '17 at 2:05
• I guess something like this: funs = Array[C[#][t] &, 4]; sols = NDSolveValue[Thread[D[funs, t] == {1, 0.72, 1.53, 1.27} funs (1 - {{1, 1.09, 1.52, 0}, {0, 1, 0.44, 1.36}, {2.33, 0, 1, 0.47}, {1.21, 0.51, 0.35, 1}}.funs)] ~Join~ Thread[(funs /. t -> 0) == ConstantArray[1, 4]], Array[C, 4], {t, 0, 1*^4}]; ParametricPlot3D[Evaluate[{sols[[1]][t], sols[[2]][t], sols[[3]][t]}], {t, 100, 1*^4}, AxesLabel -> {"x1", "x2", "x3"}, ColorFunction -> (ColorData["SunsetColors", #3] &), PlotPoints -> 200]. It seems the one who produced the Wiki image ignored the fourth coordinate for coloring. – J. M. will be back soon Apr 18 '17 at 2:44
• @J.M. Thanks, it's always interesting to see how other people do things! As for the coloring, it seems that x3 and x4 are highly correlated, which explains why you got a similar looking figure by coloring according to x3. I'm a bit curious how to color according to another function, so might post that as a separate question myself. – Chris K Apr 18 '17 at 13:54
• @Chris K Hi, sorry for bothering you again. With the help yesterday I was able to plot the biomass vs. biodensity. As I am moving to the N-species Lotka-Voltarra, I need to assign random values to $r[i]=$, so instead of saying $r[1]=RandomReal[{0, 1}]$, $r[2]=RandomReal[{0, 1}]$, $r[3]=RandomReal[{0, 1}]$ and so on, is there a commend I can use so that it assigns values to all $r[i]$? I tried using table and array, but my knowledge was too shallow. Also, if $amat$ is a matrix, what does the line $Do[a[i, j] = amat[[i, j]], (i, n), (j, n)]$ do? Thank you so much again! – Vanya Apr 18 '17 at 17:15
• @ChrisK Maybe I can sue $r[i] = RandomReal[{0, 10}, {n, 1}] // MatrixForm$ and correspond each row with an i, but how can I do that? – Vanya Apr 18 '17 at 17:19