Stability of Nonlinear Systems

Lotka-Volterra (LV) model definition

lvMod <- function(t, state, parameters){
  with(as.list(c(state, parameters)),{  
    dx <- pxx*x*x + pxy*x*y + px*x
    dy <- pyx*y*x + pyy*y*y + py*y
    list(c(dx, dy))
  })
}

LV model parameters

parameters <- c(pxx = 0.02,
                pxy = -0.23,
                px  = 0.84,
                pyx = 1.06,
                pyy = -0.04,
                py = -0.82
                )

Definition of time and initial conditions

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

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 = lvMod, parms = parameters, positive = T)
nSS$y

## x y 
## 0 0

Steady-states can also be found using simulations by allowing for a long amount of time to pass. The rootSolve package also has a function that uses the simulation approach to finding steady states. It is possible for the two approaches to identify different steady state values using the same starting point.

rSS <- runsteady(y = state, fun = lvMod, parms = parameters, times = c(0, 1e5))
rSS$y

##            x            y 
## -4.20000e+01 -4.63796e-21

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, lvMod, 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.84  0.00
## [2,] 0.00 -0.82

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 = lvMod, #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.02 and a determinant of -0.6888. The eigenvalues of this system are 0.84 and -0.82, and based on this information this steady state of this linear system is unstable and the system does not osscilate.