Stability of S-systems

S-system model definition

sMod <- function(t, state, parameters){
  with(as.list(c(state, parameters)),{  
    dx <- ax*x^gxx*y^gxy - bx*x^hxx*y^hxy
    dy <- ay*x^gyx*y^gyy - by*x^hyx*y^hyy
    list(c(dx, dy))
  })
}

S-system model parameters

parameters <- c(ax  = 1.03,
                bx  = 1,
                gxx = 0.4,
                gxy = -0.15,
                hxx = 0.2,
                hxy = 0,
                ay = 1,
                by = 1,
                gyx = 1,
                gyy = 0,
                hyx = 0,
                hyy = 1
                )

Definition of time and initial conditions

time <- seq(0, 20, 0.01)
state <- c(x = 2, #In this example the initial conditions are also used as an initial guess for the SS solver
           y = 2)

Newton-Raphson method for root solving

The package rootSolve has utilities for identifiying steady-states of ODE models defined for the deSolve package. One standard approach for finding roots of systems of equations is the Newton-Raphson method which is used here in the stode function.

nSS <- stode(y = state, fun = sMod, parms = parameters, positive = T)
nSS$y

##         x         y 
## 0.5536757 0.5536757

Linearization of nonlinear models.

To perform SS analysis of nonlinear systems, after identifying SS values, the next step is to linearize the system at the identified steady states before performing eigenvalue analysis.

A <- jacobian.full(nSS$y, sMod, parms = parameters)  #Approximation of Jacobian at SS
eigenA <- eigen(A) #Computes eigenvalues and eigenvectors of matrix A
round(A,digits = 2)

##         x     y
## [1,] 0.32 -0.24
## [2,] 1.00 -1.00

Numeric simulation to investigate stability of steady states

out <- ode(y = nSS$y*1.1 + 1,#Initial condition is based on a small perturbation from the found SS
           times = time, #Time definition
           func = sMod, #Name of function defining ODEs
           parms = parameters #Parameter values
           )

plot(out[,"time"],out[,"x"],col=4,lwd=2,type="l",xlab="Time",ylab="Values",ylim=range(c(out[,"x"],out[,"y"]))) #Plot x vs time
lines(out[,"time"],out[,"y"],col=2,lwd=2) #Plot y vs time
legend("topright",legend = c("x","y"),col = c(4,2),lwd = 2) #Define plot legend

Trace-determinant plot

traceA <- sum(diag(A))
detA <- det(A)
desA <- traceA^2 - 4*det(A)

mt <- abs(traceA)
plot(seq(-3*mt,3*mt,0.01),seq(-3*mt,3*mt,0.01)^2/4,xlab="tr A",ylab="det A",
     type="l",lwd=2,ylim=c(min(detA-10,0),max((3*mt)^2/4,detA+10)),cex.lab=1.25,cex.axis=1.25,axes=F)
axis(1, pos=0)
axis(2, pos=0)  
points(traceA,detA,pch=19,col=2,cex=1.2)

Phase-plane plot

plot(out[,"x"],out[,"y"],col=1,type="l",lwd=2,xlim=range(out[,"x"]),ylim=range(out[,"y"]),xlab="x",ylab = "y")

Stability analysis

The above defined system has a trace of -0.6790587 and a determinant of -0.0802353. The eigenvalues of this system are -0.78 and 0.1, and based on this information this steady state of this system is unstable and the system does not osscilate.