I am trying to integrate a fairly simple set of linear differential equations. When getting closer to 0, one of the variables in question (x1) starts oscillating going from negative to positive and back. This must be wrong: it’s provable that with our initial conditions it can never go below zero. I was wondering if this is an issue with the numerical algorithm for integration. The abnormal oscillations are of teh order 1e-7: much smaller than the scale of my variable.
Here’s the set of equations:
f ← function(x,p,t) {
dxdt1 = p[1]x[1](1-p[11]*x[2]-p[15]*x[4]) - p[2]*x[1]
dxdt2 = p[9]*p[2]*x[1] + p[3]x[2](1-p[12]*x[3]) - p[5]*x[2]
dxdt3 = p[5]*x[2] + p[7]x[3](1-p[13]*x[3]) - p[8]*x[3]
dxdt4 = p[10]*p[2]*x[1] + p[4]x[4](1-p[14]*x[4]) - p[6]*x[4]
return(c(dxdt1,dxdt2,dxdt3,dxdt4))
}
And I used the diffeqr package for R:
x0 = c(35.0,8.0,6.0,3.67)
sol = diffeqr::ode.solve(f, x0, tspan, p=p, saveat=0:100)
