How to solve nested ODE equations

Question

We can use the deSolve package in R for ordinary differential equations (ODE), however, I can't find a way to solve two nested ODE equations, suppose `

b'(t) = beta - k*b(t); 
a'(t) = alpha -b(t)*gamma;

where ' means differentiation. How can we solve a and b then? as a' is a function of b, we have to solve a and b simultaneously.

I got an error:

Error in lsoda(y, times, func, parms, ...) : The used combination of solvers cannot be nested.

When I tried to add the solve for b inside the ode solution for a.

Ben Bolker Ben Bolker · Accepted Answer · 2011-04-12T14:28:34

I may be confused, but you seem to be describing coupled equations, which lsoda can handle perfectly well, as follows (I implemented your ODEs but made up some parameters since I didn't know what you had in mind.)

gfun <- function(t,y,parms,...) {
  ## 'with' trick lets us write gradient in terms of variable/parameter names
  with(as.list(c(y,parms)),
       list(c(b=beta-k*b,a=alpha-b*gamma),NULL))
}

library(deSolve)
L1 <- lsoda(y=c(b=1,a=1),
            times=seq(0,10,by=0.1),
            func=gfun,
            parms=c(alpha=0.1,beta=0.2,gamma=0.05,k=0.01))

matplot(L1[,1],L1[,-1],type="l",lty=1,bty="l",las=1)

PS: this seems to be a set of coupled linear ODEs, so you should actually be able to get a full closed-form solution rather than solving them numerically. (I'm too lazy to do that right now; b(t) can be solved immediately (an "affine" equation), a(t) can be solved by integration.)

How to solve nested ODE equations

1 Answers