tests/fhn.R

In pomp: Statistical Inference for Partially Observed Markov Processes

library(pomp)

pdf.options(useDingbats=FALSE)
pdf(file="fhn.pdf")

## autonomous case
fhn <- pomp(
            data=data.frame(time=seq(0,60,by=0.1),Vobs=NA),
            times="time",
            t0=-0.1,
            skeleton.type="vectorfield",
            skeleton=function(x,t,params,...) {
              with(
                   as.list(c(x,params)),
                   c(
                     V=c*(V-V^3/3-R+i),
                     R=(V+a-b*R)/c
                     )
                   )
            }
            )

x <- array(c(0,1,1,2,1,1,0,-1),dim=c(2,2,2),dimnames=list(c("V","R"),NULL,NULL))
params <- rbind(a=c(0.7,0.5),b=c(0.8,0.5),c=c(2,5),V.0=c(1,1),R.0=c(0,0),i=c(0.8,0))
skeleton(fhn,x,t=c(0,3),params=params)
y <- trajectory(fhn,params=params,hmax=0.1)
invisible(y[,,599:601])
matplot(time(fhn),t(y["V",,]),type='l',lty=1)
plot(y[1,,],y[2,,],type='n')
points(y[1,1,],y[2,1,],pch='.',cex=3,col='black')
points(y[1,2,],y[2,2,],pch='.',cex=3,col='red')

## nonautonomous case
fhn <- pomp(
            data=data.frame(time=seq(0,20,by=0.1),Vobs=NA),
            times="time",
            t0=0,
            tcovar=seq(0,21,by=0.1),
            covar=cbind(i=sin(2*pi*seq(0,21,by=0.1))),
            skeleton.type="vectorfield",
            skeleton=function(x,t,params,covars,...) {
              c(
                V=unname(params['c']*(x['V']-(x['V']^3)/3-x['R']+covars['i'])),
                R=unname((x['V']+params['a']-params['b']*x['R'])/params['c'])
                )
            }
            )

invisible(skeleton(fhn,x,t=c(0,3),params=params))
y <- trajectory(fhn,params=params,hmax=0.01)
invisible(y[,,199:201])
matplot(time(fhn),t(y["V",,]),type='l',lty=1)
plot(y[1,,],y[2,,],type='n')
points(y[1,1,],y[2,1,],pch='.',cex=3,col='black')
points(y[1,2,],y[2,2,],pch='.',cex=3,col='red')

dev.off()