3-dimensional Plotting

Question

I am trying to visualize 3-dimensional trajectories for a system of non-linear ordinary differential equations. Building on some very valuable help I received on this forum, I have the following code:

Clear[f, g, h, p, r, l, jac, u1, u2, u3, u4, G, x, y, z, sol, xinit, \
yinit, zinit]
r = 1;    (*Recombination parameter*)
G = {{3, 3, 3, 1}, {2.5, 2.5, 
   2.5, 0.5}, {2.5, 2.5, 2.5, 0.5}, {3.5, 3.5, 3.5, 1.5}};
u1[x_, y_, z_] = 
  G[[1, 1]]*x +  G[[1, 2]]*y + G[[1, 3]]*z +  
   G[[1, 4]]*(1 - x - y - z) ;
u2[x_, y_, z_] = 
  G[[2, 1]]*x +  G[[2, 2]]*y + G[[2, 3]]*z +  
   G[[2, 4]]*(1 - x - y - z);
u3[x_, y_, z_] = 
  G[[3, 1]]*x +  G[[3, 2]]*y + G[[3, 3]]*z +  
   G[[3, 4]]*(1 - x - y - z) ;
u4[x_, y_, z_] = 
  G[[4, 1]]*x +  G[[4, 2]]*y + G[[4, 3]]*z +  
   G[[4, 4]]*(1 - x - y - z) ;
ualpha[x_, y_, 
   z_] = (x*u1[x, y, z]) + (y*u2[x, y, z]) + (z*
     u3[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);
us[x_, y_, z_] = (x*u1[x, y, z]) + (y*u2[x, y, z]);
ua[x_, y_, z_] = (z*u3[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);
uc[x_, y_, z_] = (x*u1[x, y, z]) + (z*u3[x, y, z]);
ud[x_, y_, z_] = (y*u2[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);
F1[x_, y_, 
   z_] = ((1 - r)*x*u1[x, y, z]/ualpha[x, y, z]) + (r*us[x, y, z]*
     uc[x, y, z]/((ualpha[x, y, z])^2)) - x;
F2[x_, y_, 
   z_] = ((1 - r)*y*u2[x, y, z]/ualpha[x, y, z]) + (r*us[x, y, z]*
     ud[x, y, z]/((ualpha[x, y, z])^2)) - y;
F3[x_, y_, 
   z_] = ((1 - r)*z*u3[x, y, z]/ualpha[x, y, z]) + (r*ua[x, y, z]*
     uc[x, y, z]/((ualpha[x, y, z])^2)) - z;
nmax = 100;
tmax = 100;
func1 = {};
func2 = {};
For[k = 1, k <= nmax, k++, num1 = RandomReal[];
 num2 = RandomReal[];
 max = Max[num1, num2];
 min = Min[num1, num2];
 rand1 = RandomReal[{0, min}];
 rand2 = RandomReal[{min, max}];
 rand3 = RandomReal[{max, 1}];
 solution = 
  NDSolve[{x'[t] == F1[x[t], y[t], z[t]], 
    y'[t] == F2[x[t], y[t], z[t]], z'[t] == F3[x[t], y[t], z[t]], 
    x[0] == rand1, y[0] == rand2 - rand1, z[0] == rand3 - rand2}, {x, 
    y, z}, {t, 0, tmax}];
 plotfunc1 = 
  ParametricPlot3D[{x[t], y[t], z[t]} /. solution, {t, 0, tmax}, 
    PlotRange -> All, BaseStyle -> Arrowheads[{0, .015, 0.015, 0}], 
    AxesLabel -> {"x", "y", "z"}] /. Line -> Arrow;
 plotfunc2 = 
  ParametricPlot3D[{x[t], y[t], z[t]} /. solution, {t, 0, tmax}, 
   PlotRange -> All, ColorFunction -> Hue, PlotTheme -> "Marketing", 
   AxesLabel -> {"x", "y", "z"}];
 AppendTo[func1, plotfunc1];
 AppendTo[func2, plotfunc2]]
Show[func1]
Show[func2]

The code does the following. It randomly picks 100 points (as initial conditions) from the $xyz$ space such that $x,y,z\geq 0$ and $x+y+z \leq 1$ and plots the trajectories for the system: $\dot{x} = F_1(x,y,z),$ $\dot{y} = F_2(x,y,z)$ and $\dot{z} = F_3(x,y,z).$ I used to two different plotting styles whose output are attached. As one can see in the figures there is a lot of unused space which am unable to get rid of. Any help on removing the unused space and improving the clarity of figures will be greatly appreciated.

Answer 1

The code does the following. It randomly picks 100 points (as initial conditions) from the $xyz$ space such that $x,y,z \geq 0$ and $x+y+z \leq 1$ and plots the trajectories for the system: $\dot{x} = > F_1(x,y,z),$ $\dot{y} = F_2(x,y,z)$ and $\dot{z} = F_3(x,y,z).$

I would ideally want to have a bunch of random initial condition points on the pyramid and then trace their trajectory.

Your code doesn't appear to randomly sample your space and definitely doesn't ensure your points start on the pyramid face ($x+y+z=1$). Unsure if that's important for you. Either way, you should try to avoid loops in Mathematica - opting for Map, Table, ...etc. when possible. I've rewritten your code to be a little more "Mathematica-y". It also processes everything 10x faster on my machine when nmax=100. That doesn't really matter here since the plotting takes way more time than the processing - but best to do things right because they might matter next time!

As for the pyramid shaped plot - I've thrown together a quick solution. Basically you just remove the bounding box and replace it with transparent faces to define your pyramid. You could do a lot to clean it up further, but hopefully this gives you what you need.

Clear[vec, F, u, v, p, sol];
vec[p_] := Flatten[{p, 1 - Total[p]}];

nmax = 100;
tmax = 100;
r = 1;(*Recombination parameter*)
G = {{3, 3, 3, 1}, {2.5, 2.5, 2.5, 0.5}, {2.5, 2.5, 2.5, 0.5}, {3.5, 3.5, 3.5, 1.5}};

i1 = {All, {1, 2}, {3, 4}, {1, 3}, {2, 4}}; (* indices for dot product giving ualpha - ud *)
i2 = {{2, 4}, {2, 5}, {3, 4}};(* indices for products like us*uc *)

u[p_] := (#.vec[p] & /@ G); (* u1 - u4 *)

v[p_] := Total[((vec[p] u[p])[[#]])] & /@ i1; (* ualpha - ud *)

p[p_] := r Table[v[p][[i2[[k, 1]]]]*v[p][[i2[[k, 2]]]], {k, 3}]/(v[p][[1]]^2);(* 2nd term in F1-F3 *) 

F[x_, y_, z_] = With[{w = {x, y, z}}, (1 - r)*w . Rest[u[w]]/v[w][[1]] + p[w] - w]; (* F1 - F3 *)

region = ImplicitRegion[x + y + z <= 1 && x >= 0 && y >= 0 && z >= 0, {x, y, z}]; (* better way to get uniform points *)
{xs, ys, zs} = Thread[#/Total[#] & /@ RandomPoint[region, nmax]] (* divide by total so every point starts on pyramid *);

solution = First@NDSolve[{x'[t] == F[x[t], y[t], z[t]][[1]],y'[t] == F[x[t], y[t], z[t]][[2]],z'[t] == F[x[t], y[t], z[t]][[3]], x[0] == xs, y[0] == ys, z[0] == zs}, {x, y, z}, {t, 0, tmax}];

sol[t_] := Thread[Through[({x, y, z} /. solution)[t]]];

p1 = ParametricPlot3D[sol[t], {t, 0, 100},PlotRange -> {{0, 1}, {0, 1}, {0, 1}},BaseStyle -> Arrowheads[{0, .015, 0.015, 0}], 
  ColorFunction -> "Rainbow", Boxed -> False, AxesLabel -> {x, y, z},AxesStyle -> Directive[Black, Bold, Thick, 16], 
  AxesEdge -> {{-1, -1}, {-1, -1}, {-1, -1}}, Ticks -> ConstantArray[Table[k/4.,{k, 4}],3]];(* plot without box and axes along pyramid *)
p2 = Graphics3D[{Black, Opacity[0.1], HalfPlane[{{1, 0, 0}, {0, 1, 0}}, {0, -1, 0}],HalfPlane[{{0, 0, 1}, {0, 1, 0}}, {0, -1, 0}], HalfPlane[{{0, 0, 1}, {1, 0, 0}}, {-1, 0, 0}]}]; (* pyramid faces *)
Show[p1, p2]