Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}]
share|improve this question
    
First, use := and patterns to define functions. Second, without the function m, I doubt someone will understand what the issue is. –  Michael E2 Jul 19 '13 at 0:32
add comment

2 Answers

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 (say it approaches asymptotically for example). Don't despair though; you can still get numerical solutions within NDSolve and, say, track equilibrium quantitatively, incorporating an equilibrium threshold in the form of a Method within the numerical integrator NDSolve that then Reaps the 'synthetic data' stating just how small the derivatives of your state variables (x, y, z here) 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, with true scientific sentiment, 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 even tells it to automatically stop the integration automatically once the relative equilibria time is found.

I assure you this can be done as I have done it with 100+ coupled differential equations within a Manipulate that isn't too slow on my computer. I must admit to having done exactly what you are now: thinking I could input the equilibrium condition as zeroes and then just tell the machine to Solve. A few nine to twelve hour waits while Mathematica populates all 8 GB of ram and then erred informed me of my folly...

If you want to begin to format your equations for input into NDSolve, the answer to this may be of interest to you. In fact, while this was a novice question (& duplicate), in the question, when posting it, I provided the code to constructing the Manipulate to vary the physical constants etc. via real-time sliders as I mentioned earlier, which more than encompasses the proper formatting necessary for input into NDSolve and etc., after which 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. I'm more than willing to try to help you overcome any obstacles (perhaps in chat) if you do think this seems a profitable route. Sorry I couldn't be more detailed here, the surprise of unexpectedly proposed grad school in only a few weeks has me in a tither. If I think of further details I'll update the answer.

share|improve this answer
add comment

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 ?

share|improve this answer
    
First I think that it's Clear and not clear. –  Öskå Jul 19 '13 at 8:40
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.