# Linear Stability Analysis of dynamical system

I'm starting to use Mathematica for some linear stability analysis of a discrete non-linear dynamical system. I can't find an on-line tutorial for it, and I'm quite at a loss (in fact I can't even replicate what I've already done on my paper notebook...). Here's my progress so far:

1. I define the function I am studying (2 variables, 2 dimensions). This works.
2. Calculate the Jacobian Matrix in a generic point. This works.
3. Calculate eigenvalues in a generic point. This works.
4. I want to define a rule f[x,y] -> {x,y} and apply it to the Jacobian before calculating eigenvalues. This doesn't work. I've tried different variations of this, even writing the actual expression, and the substitution is pretty obvious (the function has a bunch of exponentials so the derivative is close to the original function). So what I do is Jac /. f[x,y] -> {x,y} (or variations where I actually write the expression) but it doesn't seem to do anything.
5. Once I had these eigenvalues at the solution points, I want to see if their absolute value is bounded by 1 given a condition on the parameters. I have no idea how to show it. The system depends on 6 parameters and the condition for stability is non-linear and non-nice. In this case I know the formula so my question is how to tell Mathematica "this is positive, this is less than that, etc." to see if the eigenvalue is greater or smaller than 1.

I can't solve f[x,y] == {x,y} (or I don't know how to do it at least), and there are two points that I can check manually ({0,0} and {0,k}) but there is a third positive solution that can't be found analytically. That's why I need to use that rule and not simply substitute x and y by the solutions.

Bonus questions (i.e. I don't need the answer but if you could point me to a manual or something I'd appreciate that): I can prove by hand that there is a single other solution, but how would you go about it in Mathematica? Also, is there a way to go the other way around in the eigenvalues? I mean, tell Mathematica "The eigenvalues are bounded by 1, what conditions does this impose on parameters?".

Edit: Here's the actual code. Doubt it will be very useful.

g[H_, P_] := {Subscript[s, p]*P + (1 - E^(-P))*H*d,
k/(k + H)*(Subscript[s, h] + f*E^(-P) + f*(1 - E^(-P))*(1 - d))*H}
Jac := D[g[H, P], {{H, P}}] // MatrixForm
Jac /. H*k*(E^(-P)*f + (1 - d)*(1 - E^(-P))*f + Subscript[s, h])/(H +
k) -> H


Last step doesn't do anything. Neither any variation of this.

Edit2: I have removed MatrixForm from the code. Same thing happens. I also removed subscripts just in case. Thank you for your help!

• We cannot say anything without the actual code that you have used... – Henrik Schumacher Feb 28 '18 at 21:15
• @HenrikSchumacher added code. – bernatguillen Feb 28 '18 at 21:20
• do not use subscripts. Not good for coding. Also not good idea to attach MatrixForm to a result. Since now you can't work with the result any note since MatrixForm is wrapper around it. That is why nothing happens in last step. – Nasser Feb 28 '18 at 21:32
• @Nasser thanks! I removed MatrixForm but the same happens. It does work if I do P-> 0 /. H->0 (which is one of the solutions). – bernatguillen Feb 28 '18 at 21:36
• What's the idea behind step 4 -- what is f[x,y]? In your code, what's the purpose of the last line? Could you give a little background on the model? – Chris K Mar 1 '18 at 0:01