4
$\begingroup$

Sol is the evaluation at $t=100$ of the solution of a nonlinear system of three differential equations with initial values at $a$, $b$, $c$. This is obtained with the function ParametricNDSolveValue. 

 sol = ParametricNDSolveValue[ {x'[t] == -y[t], y'[t] == -z[t], z'[t] == -0.915 x[t]+(1-0.915 x[t]^2) y[t]-z[t], x[0] == a, y[0] == b, z[0] == c}, {x[100], y[100], z[100]}, {t, 0, 100}, {a, b, c} ]

The point $\{xeq,yeq,zeq\}=\{0,0,0\}$ is one of the stable equilibriums of the alleged system. I would like to obtain the stable region of the this equilibrium so for that I'm trying the following: calculate the distance from a point evaluated with Sol to the equilibrium and then plot all these points. The points I evaluate is the points of the unit cube obtained with a step of 0.1 in each of the coordinates:

Table[EuclideanDistance[sol[a, b, c], {xeq, yeq, zeq}], {a,0, 1,0.1}, {b, 0, 1, 0.1},{c,0,1,0.1}]

This provides me with the distance. The problems are two:

1) I don't know how to obtain back the point $\{a,b,c\}$ (the argument of Sol) for the distances that are less than, say, 0.1 from the equilibrium.

2) I actually have only tried with only one or two variables in a small range, say,

Table[EuclideanDistance[sol[a, b, 0.1], {xeq, yeq, zeq}], {a,0.4, 0.6,0.1}, {b, 0.7, 0.9, 0.1}]

because I usually get the following 

General: 1/4.60387*10^307 is too small to represent as a normalized machine number; precision may be lost.

which makes the evaluation last forever not showing any outcome.

How can I overcome these issues?

Improve this question
$\endgroup$
7
  • 3
    $\begingroup$ You need to include the defintiion of sol in your question. $\endgroup$ Commented Jun 14, 2019 at 8:27
  • $\begingroup$ @KraZug, ok I included it. $\endgroup$ Commented Jun 14, 2019 at 8:46
  • $\begingroup$ And {xeq, yeq, zeq}? $\endgroup$ Commented Jun 14, 2019 at 12:23
  • $\begingroup$ Your equations have not converged to a steady state at $t=100$ for $a=1,b=1,c=1$. In fact, you have $z \to infty$ as $t \to \infty$. Perhaps you need to check the equation you have for $z$. $\endgroup$ Commented Jun 14, 2019 at 12:29
  • $\begingroup$ @KraZug, yes, you're right, sorry I had several systems saved as labels. I corrected the system. $\endgroup$ Commented Jun 14, 2019 at 18:36

2 Answers 2

Reset to default
4
$\begingroup$

Here's a way rather in line with the OP's approach:

(* returns the squared distance to the origin *) sol = ParametricNDSolveValue[{x'[t] == -y[t], y'[t] == -z[t], z'[t] == -0.915 x[t] + (1 - 0.915 x[t]^2) y[t] - z[t], x[0] == a, y[0] == b, z[0] == c, WhenEvent[ x[t]^2 + y[t]^2 + z[t]^2 < Min[(a^2 + b^2 + c^2)/100, 10^-6], "StopIntegration"]}, Indexed[x["ValuesOnGrid"], -1]^2 + (* This is the same idea as KraZug's *) Indexed[y["ValuesOnGrid"], -1]^2 + (* use of the last value, but with a *) Indexed[z["ValuesOnGrid"], -1]^2, (* different code *) {t, 0, 100}, {a, b, c}]; stab = RegionPlot3D[ Quiet@sol[a, b, c] < 1 (* a smaller threshold is not really necessary *) {a, -2, 2}, {b, -2, 2}, {c, -2, 2}, Mesh -> None]; // AbsoluteTiming (* {29.8274, Null} <-- Takes a while! *)

We'll show the region with the flow indicated on the surface, which should be tangent to it (to the gold part, not the blue part). We can see that the stability region extends further out the top and bottom:

(* lazy way to get the vector field *) {state} = NDSolve`ProcessEquations[{x'[t] == -y[t], y'[t] == -z[t], z'[t] == -0.915 x[t] + (1 - 0.915 x[t]^2) y[t] - z[t], x[0] == 1, y[0] == 2, z[0] == 3}, (* random IC, not used *) {x, y, z}, {t, 0, 100} ]; nf = state@"NumericalFunction"; Show[Cases[stab, GraphicsComplex[p_, ___] :> VectorPlot3D[nf[0, {a, b, c}], {a, -2, 2}, {b, -2, 2}, {c, -2, 2}, VectorPoints -> p, VectorScale -> {0.045, Automatic, None}], Infinity], stab /. RGBColor[r_, g_, b_] :> RGBColor[r, g, b, 0.6] (* opacity was an afterthought *) ]

enter image description here

If you want to use a greater number of PlotPoints in RegionPlot, you can downsample the points for VectorPlot with something like the following:

VectorPoints -> Union[p, SameTest -> (EuclideanDistance[##] < 0.2 &)]
Improve this answer
$\endgroup$
4
$\begingroup$

So the issue is that at some parameter the solution diverges, which then gives vectors that are too large to use EuclideanDistance on. So we need to catch that, for instance in this way:

sol = ParametricNDSolveValue[{x'[t] == -y[t], y'[t] == -z[t], z'[t] == -0.915 x[t] + (1 - 0.915 x[t]^2) y[t] - z[t], x[0] == a, y[0] == b, z[0] == c}, {x, y, z}, {t, 0, 10000}, {a, b, c}] tab = Flatten[Table[cursol = sol[a, b, c]; endpt = cursol[[1, 1, 1, 2]]; {a, b, c, EuclideanDistance[Through[cursol[endpt]], {0, 0, 0}]}, {a, 0, 1, 0.1}, {b, 0, 1, 0.1}, {c, 0, 1, 0.1}], 2];

Here sol now returns the whole interpolation function, and cursol[[1, 1, 1, 2]] digs into the interpolation functions to return the endpoint of the integration. The table returns a list {x,y,z,f}, which we can split based on the values of f:

ListPointPlot3D[{Select[tab, #[[4]] < 0.01 &][[All, 1 ;; 3]], Select[tab, #[[4]] >= 0.01 &][[All, 1 ;; 3]]}, PlotStyle -> {Black, Red}]

ListPointPlot3D with black points representing converged solutions

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.