Autonomous equations' critical points and eigenvalues

Question

I have three-dimensional autonomous equations and want to find their critical points and corresponding eigenvalues. They have 6 critical points each.

x, y and z are functions of time and m is a function of z/y

x'=-1 - z - 3*y + x^2 - x*z,
y'=(x*z/m[z/y]) - y*(2*z - 4 - x),
z'=-(x*z/m[z/y]) - 2 z*(z - 2),

I wrote this code but it didn't work. Top of the page says running but nothing happens.

f[x, y, z] == -1 - z - 3*y + x^2 - x*z
g[x, y, z] == (x*z/m[z/y]) - y*(2*z - 4 - x)
h[x, y, z] == -(x*z/m[z/y]) - 2 z*(z - 2)

sys == {D[x[t], t] == f[x[t], y[t], z[t]], 
  D[y[t], t] == g[x[t], y[t], z[t]], 
  D[z[t], t] == h[x[t], y[t], z[t]]};

Solve[{f[x, y, z] == 0, g[x, y, z] == 0, h[x, y, z] == 0}, {x, y, z}]

Answer 1

For me, Solve reports an error until only the right side is extracted (without the derivative) via

eqns =
 {x' == -1 - z - 3*y + x^2 - x*z,
  y' == (x*z/m[z/y]) - y*(2*z - 4 - x),
  z' == -(x*z/m[z/y]) - 2 z*(z - 2)}

and isolating RHS to avoid a "This cannot be solved via the methods available toSolve" error,

listRHS = eqns /. (lhs_ == rhs_) :> (rhs)    

Solve[Thread[listRHS == 0], {x, y, z}]

which together make the seperate definition of f[x,y,z], etc. unnecessary for a moment (which may end up actually being more convenient for you given that Solve runs seemingly indefinitely).

The reason for this is because some systems never reach true equilibrium (asymptotically for example). You can still get numerical solutions within NDSolve and track equilibrium, incorporating a threshold in the form of a Method within the numerical integrator that then Reaps the just how small the derivatives of your state variables are at that time.

For example if

dstateVariables = {x'[t], y'[t], z'[t]}

then your addition to Sjoerd's answer here (while they are all excellent, his lack of 3rd party functions supports things like PrecisionalGoal, AccuracyGoal etc.) would be:

equilibriaThreshold = 
 10^-11.1;  (* Units of state variable input per unit time. *)
equiCond = Total@Thread[Abs@dstateVariablesExplicit];

Method -> {"EventLocator", "Event" -> equiCond < equilibriaThreshold, 
"EventAction" :> Throw[end = t, "StopIntegration"]}

Where as you can see, the equilibrium condition achieved in my system (which is the sum of 100+ state variables' derivatives) is damned small , so pretty much zero (steady-state). The Throw term in the last line tells it to automatically stop the integration automatically once the relative equilibria time is found.

This works with 100+ coupled differential equations within a Manipulate that isn't too slow. Similarly, I tried using Solve with the steady-state scenario; a few nine to twelve hour waits while Mathematica populates all 8 GB of RAM were my punishment then.

If you want to begin to format your equations for input into NDSolve, the answer to this may be of interest. In fact, while this was a novice question (& duplicate), in the question, when posting it, I provided the code (mostly found elsewhere) to construct the Manipulate to vary any physical constants etc. via real-time sliders, which more than encompasses the proper formatting necessary for input into NDSolve etc. Then, all you've left to do is step through the motions of other answers (such as Sjoerd's linked earlier) to come up with the equilibria.

I just started about four months ago, so it's certainly within your capacity to set this up. If I think of further details I'll update the answer.

Answer 2

I wrote it in another way and got an answer, but under the assumption that m = constant :

Clear[x, y, z];
crit = Solve[
  {-1 - z - 3 y + x^2 - x*z == 0, 
   (x*z)/m + 4 y - 2 y*z + x*y == 0,
   -(x*z)/m - 2 z^2 + 4 z == 0},
  {x, y, z}]

{{x -> -4, y -> 5, z -> 0}, {x -> -1, y -> 0, z -> 0}, 
  {x -> 0, y -> -1, z -> 2}, {x -> 1, y -> 0, z -> 0}, 
  {x -> (3 m)/(1 + m), y -> 1/6 (-6 + (9 m)/(1 + m)^2 + 
      (18 m^2)/(1 + m)^2 + 3/(1 + m) - (12 m)/(1 + m)), 
   z -> (1 + 4 m)/(2 (1 + m))}, 
  {x -> -((2 (-1 + m))/(1 + 2 m)), 
   y -> -((6 m - (4 (-1 + m)^2)/(1 + 2 m)^2 - (8 (-1 + m)^2 m)/(1 + 2 m)^2 +
      (2 (-1 + m))/(1 + 2 m) - (8 (-1 + m) m)/(1 + 2 m))/(6 m)), 
   z -> (-1 + 3 m + 4 m^2)/(m (1 + 2 m))}}

but if take m(z/y) after some time Mathematica returns:

No more memory available.
Mathematica kernel has shut down.
Try quitting other applications and then retry.

What do you think guys ?