# R/euler-RK4-function.R In seem: Simulation of Ecological and Environmental Models

```euler <- function(x0, t, f, p, dt){

# Euler numerical solution of model
# arguments: x0 initial condition
#            t times for output
#            f model
#            p parameters
#            dt time step
# arrays to store results
nt <- length(t); nx <- length(x0)
ns <- floor((t[2]-t[1])/dt)
X <- matrix(nrow=nt,ncol=nx)
# first value of X is initial condition
Xt <- x0; tt <- t[1]
# loops
for(i in 1:nt){
X[i,] <- Xt
for(j in 1:ns){
tt <- tt + dt;Xt <- Xt + dt*f(tt,p,Xt)
# X assumed positive
for(k in 1:nx) if(Xt[k]<0) Xt[k] <- 0
}
} # end of run
return(X)
} # end of function

RK4 <- function(x0, t, f, p, dt){

# RK4 numerical solution of model
# arguments: x0 initial condition
#            t times for output
#            f model
#            p parameters
#            dt time step
# arrays to store results
nt <- length(t); nx <- length(x0)
ns <- floor((t[2]-t[1])/dt)
X <- matrix(nrow=nt,ncol=nx)
# first value of X is initial condition
Xt <- x0; tt <- t[1]
# loops
for(i in 1:nt){
X[i,] <- Xt
for(j in 1:ns){
k1<- f(tt,p,Xt)
k2<- f(tt+dt/2,p,(Xt+k1*dt/2))
k3<- f(tt+dt/2,p,(Xt+k2*dt/2))
k4<- f(tt+dt,p,(Xt+k3*dt))
kavg <- (dt/6)*(k1+ 2*k2+ 2*k3+ k4)
tt <- tt + dt; Xt <- Xt + kavg
# X assumed positive
for(k in 1:nx) if(Xt[k]<0) Xt[k] <- 0
}
} # end of run
return(X)
} # end of function
```

## Try the seem package in your browser

Any scripts or data that you put into this service are public.

seem documentation built on April 14, 2017, 9:12 p.m.