1
$\begingroup$

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.Figure 2

Figure 1

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Have you tried changing the PlotRange option in ParametricPlot3D? Perhaps try something like PlotRange -> {{0,1},{0,0.2},{0,1}} $\endgroup$ Commented Jun 23, 2021 at 15:56
  • $\begingroup$ @bRost03: I tried using the option you suggested but it does not cover the whole space. The whole space I intend to cover is $x,y,z \geq 0$ and $x+y+z \leq 1.$ $\endgroup$ Commented Jun 23, 2021 at 16:17
  • $\begingroup$ So you want the plot to be some sort of pyramid shape? $\endgroup$ Commented Jun 23, 2021 at 16:22
  • $\begingroup$ @bRost03: Yes, something like that (pyramid shape). I would ideally want to have a bunch of random initial condition points on the pyramid and then trace their trajectory. $\endgroup$ Commented Jun 23, 2021 at 16:27
  • $\begingroup$ @egt123, thus you should use the RegionFunctions with you formulas of volume $\endgroup$ Commented Jun 23, 2021 at 16:29

1 Answer 1

Reset to default
3
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ This gives almost everything I need. Thanks a lot! $\endgroup$ Commented Jun 23, 2021 at 21:26

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.