[?] Find the stationary states of an EDO system modelling the yellow fever?

Questions
, ,
1

I am modeling the yellow fever for an EDO system, in which I have to find the stationary states of the system, and I have great difficulty in finding something that can be very simple. I put the program to run and it keeps running forever. I thought it was the model, so I compiled another model in which I had already obtained the stationary states quickly, and the same one had the same problem. I tested on two other machines and the problem persists, including for the model I had compiled before and had got the points. I suppose it is error writing the system in wolfram to get the points of equilibrium, but the point is that this writing format I had already used to ober stationary states of the ODE in wolfram. What is wrong with calling the equations? Please tell me if you need the written template to try to help and excuse me for the English and if it is a trivial question, try to relieve me because I am new to the software.

Solve[-a*S[t]*P[t] + b*R[t] - c*S[t] + d*(S[t] + I1[t] + R[t]) == 0 &&  a*S[t]*P[t] - e*I1[t] - (I1[t]*(c + P[t])) == 0 &&  
e*I1 - b*R[t] ==  0 && -f*N1[t]*I1[t] + g*(N1[t] + P[t])*(1 - (N1[t] + P[t])/k) - h*N1[t] - c*N1[t] == 0 && 
  f*N1[t]*I1[t] - h*P[t] - c*P[t] == 0, {S[t], I1[t], R[t], N1[t], P[t]}, t]
2

There is an issue with this formulation in that the ,t at the end should not be there (DSolve needs an independent variable, Solve has no notion of that). Also I1 appears in one place without the t dependency; I assume that was a mistake. But it is a difficult system to handle with all the symbolic parameters. It might be more effective to only solve when given numeric values for {a,b,...}. Could be done as below.

exprs = {-a*S[t]*P[t] + b*R[t] - c*S[t] + d*(S[t] + I1[t] + R[t]), 
   a*S[t]*P[t] - e*I1[t] - (I1[t]*(c + P[t])), 
   e*I1[t] - b*R[t], -f*N1[t]*I1[t] + 
    g*(N1[t] + P[t])*(1 - (N1[t] + P[t])/k) - h*N1[t] - c*N1[t], 
   f*N1[t]*I1[t] - h*P[t] - c*P[t]};

params = {a, b, c, d, e, f, g, h, k};
vars = {S[t], I1[t], R[t], N1[t], P[t]};

solns[a0_?NumberQ, b0_?NumberQ, c0_?NumberQ, d0_?NumberQ, e0_?NumberQ,
   f0_?NumberQ, g0_?NumberQ, h0_?NumberQ, k0_?NumberQ] := 
 NSolve[exprs /. 
   Thread[params -> {a0, b0, c0, d0, e0, f0, g0, h0, k0}]]

Example:

 solns[.1, .3, .2, .5, .6, .4, .7, .11, .6]

(* Out[460]= {{I1[t] -> -0.833013435701, N1[t] -> -4.46571428571, 
  P[t] -> 4.8, R[t] -> -1.6660268714, 
  S[t] -> -9.71849008317}, {I1[t] -> -0.775, N1[t] -> -4.8, 
  P[t] -> 4.8, R[t] -> -1.55, 
  S[t] -> -9.04166666667}, {I1[t] -> -0.775, N1[t] -> 0.5, 
  P[t] -> -0.5, R[t] -> -1.55, 
  S[t] -> 4.65}, {I1[t] -> -0.464469178082, N1[t] -> 0.834285714286, 
  P[t] -> -0.5, R[t] -> -0.928938356164, 
  S[t] -> 2.78681506849}, {I1[t] -> 0., N1[t] -> 0.334285714286, 
  P[t] -> 0., R[t] -> 0., S[t] -> 0.}, {I1[t] -> 0., N1[t] -> 0., 
  P[t] -> 0., R[t] -> 0., S[t] -> 0.}} *)
3

Hi,

I think that there are a some issues here. I suppose that you want to set the rhs of ODEs to zero. In that case we are looking for fixed values of S,I1,R,N1 and P to solve the equations. If you knew the parameters in the equations the solution is quite straight forward and fast:

Solve[-a*S*P + b*R - c*S + d*(S + I1 + R) == 0 && a*S*P - e*I1 - (I1*(c + P)) == 0 && 
e*I1 - b*R == 0 && -f*N1*I1 + g*(N1 + P)*(1 - (N1 + P)/k) - h*N1 - c*N1 == 0 && f*N1*I1 - h*P - c*P == 0 /. {a -> 1, b -> 1, c -> 1, d -> 1, e -> 1, f -> 1, g -> 1, h -> 1, k -> 1}, {S, I1, R, N1, P}]
(*{{I1 -> 0, R -> 0, N1 -> -1, P -> 0}, {I1 -> 0, R -> 0, N1 -> 0, 
P -> 0}, {S -> -3, I1 -> -1, R -> -1, N1 -> -2, P -> 1}, {S -> -6, 
I1 -> -2, R -> -2, N1 -> -1, P -> 1}}*)

Note that I removed the time dependence of the variables. For these parameters it correctly identifies the origin as a fixed point; I suppose that the other fixed points are irrelevant in your modelling situation, because they are negative.

You can explicitly add positivity conditions like so:

Solve[(-a*S*P + b*R - c*S + d*(S + I1 + R) == 0 && a*S*P - e*I1 - (I1*(c + P)) == 0 && 
e*I1 - b*R == 0 && -f*N1*I1 + g*(N1 + P)*(1 - (N1 + P)/k) - h*N1 - c*N1 == 0 &&
f*N1*I1 - h*P - c*P == 0 && S >= 0 && I1 >= 0 && R >= 0 && N1 >= 0 && P >= 0) /. {a -> 1, b -> 1, c -> 1, d -> 1, e -> 1, 
f -> 1, g -> 1, h -> 1, k -> 1}, {S, I1, R, N1, P}]
(*{{I1 -> ConditionalExpression[0, S >= 0], 
  R -> ConditionalExpression[0, S >= 0], 
  N1 -> ConditionalExpression[0, S >= 0], 
  P -> ConditionalExpression[0, S >= 0]}}*)

In theory you could also try to use functions like Reduce to solve the system. Given the large number of parameters and nonlinear expressions I guess that the general conditions are quite complex.

Reduce[-a*S*P + b*R - c*S + d*(S + I1 + R) == 0 && a*S*P - e*I1 - (I1*(c + P)) == 0 && 
e*I1 - b*R == 0 && -f*N1*I1 + g*(N1 + P)*(1 - (N1 + P)/k) - h*N1 - c*N1 == 0 && f*N1*I1 - h*P - c*P == 0 && S > 0 && I1 > 0 && R > 0 && N1 > 0 && P > 0, {S, I1, R, N1, P}, Reals]

hasn’t finished on my machine yet. Do you have a priori values for the parameters?

Cheers,

Marco

4

Sorry, your post hasn’t shown up until I replied. The system sometimes has quite a delay… Perhaps I am sitting on the wrong side of the pond…

Cheers,

Marco