1

I have a system of three ODEs that I want to find the steady states of, and whether or not they're stable via a simple analysis. To do so, I just calculate the eigenvalues of the Jacobian evaluated at each steady state, and then determine what conditions (based on my parameters) are required to have a negative eigenvalue.

That works fine for the most part, but in my original analysis I did not include the parameter $\kappa$, equivalent to setting $\kappa=1$ in the supplied code. When I force $\kappa$ to be any particular variable, I can see how my conditions change (easily noted when $\kappa=\pi$, as $\pi$ doesn't show up in the analysis otherwise), but if I don't specify a $\kappa$ then the code tells me that it's impossible (the getConditions function returns False)

Any idea why that's happening?

Also, I'm sure there's a more proper way of doing this analysis with more built in functions, and declaring my ODEs properly, but I'm pretty new to Mathematica.

κ = π
(* Clear[κ] (* uncomment to see the problem *)
δS = F S/(S + R) - γ S P - κ S;
δR = F R/(S + R) χ - κ R;
δP = (A γ - ν) S P - ν P R - κ P;

steadyStates = 
  Assuming[0 < χ < 1 && 0 < A < 1 && γ > 0 && ν > 0 && 
    0 < F <= 1 && 1 <= κ <= 10 && A γ > ν,
   Solve[{0 == δS, 0 == δR, 0 == δP}, {S, R, P}]];
J = Simplify[({
     {\!\(
       \*SubscriptBox[\(∂\), \(S\)]δS\), \!\(
       \*SubscriptBox[\(∂\), \(R\)]δS\), \!\(
       \*SubscriptBox[\(∂\), \(P\)]δS\)},
     {\!\(
       \*SubscriptBox[\(∂\), \(S\)]δR\), \!\(
       \*SubscriptBox[\(∂\), \(R\)]δR\), \!\(
       \*SubscriptBox[\(∂\), \(P\)]δR\)},
     {\!\(
       \*SubscriptBox[\(∂\), \(S\)]δP\), \!\(
       \*SubscriptBox[\(∂\), \(R\)]δP\), \!\(
       \*SubscriptBox[\(∂\), \(P\)]δP\)}
    }), Assumptions -> 
    S >= 0 && R >= 0 && P >= 0 && 0 < χ < 1 && 
     0 < A < 1 && γ > 0 && ν > 0 && 0 < F <= 1 && 
     1 <= κ <= 10 && A γ > ν];
evs = Map[Simplify[Eigenvalues[J /. #1]] &, steadyStates];

getConditionsForStability[eigenvalue_] := 
 Simplify[Reduce[eigenvalue < 0, {F}],
  Assumptions -> 
   0 < χ < 1 && 0 < A < 1 && γ > 0 && ν > 0 && 
    0 < F <= 1 && 1 <= κ <= 10 && A γ > ν] 

getConditionsForStability /@ evs[[4]]
getConditionsForStability[evs[[4]]]
m_goldberg
  • 107,779
  • 16
  • 103
  • 257
Matt
  • 11
  • 1
  • Curious, I've rewritten the code several ways, and they all do the same thing. Could you post the actual ODE you need to solve? – Feyre Jun 06 '16 at 11:50
  • It's just the same three equations in the top, with respect to time. $\frac{dS}{dt} = F S/(S + R) - γ S P - κ S$, $\frac{dR}{dt} = F R/(S + R) χ - κ R$, $\frac{dP}{dt} = (A γ - ν) S P - ν P R - κ P$. With the parameters following the rules that I had copied and pasted a bunch of times in the code, except $\kappa$ which should be able to be $1\leq\kappa$ – Matt Jun 06 '16 at 13:31
  • 1
    My guess is that kappa cannot be determined based on the other variables. It's definitely not a code problem. Defining the ODE explicitly seems too heavy for mathematica. – Feyre Jun 06 '16 at 15:04

0 Answers0