# ======================================================== #
# Numerical Integration of Ordinary Differential Equations #
# ======================================================== #


# This library contains a didactical implementation of a numerical integrator
# of a system of ordinary differential equations.
# Note that this implementation is intended to demonstrate how such 
# techniques can be implemented in R. It does not represent the state of the 
# art of numerical integration of such differential equations.


# created and maintained by 
# Peter Reichert
# EAWAG
# Duebendorf
# Switzerland
# reichert@eawag.ch


# First version: Dec.  22, 2002
# Last revision: Jan.  07, 2007


# Overview of functions
# =====================

# ode:            numerical integration of deterministic ordinary differential 
#                 equations


# =========================================================================== #


# numerical integration of ordinary differential equations
# ========================================================

ode <- function(rhs,x.ini,param,t.out,dt.max=NA,algorithm="euler",...)
{
   # -----------------------------------------------------------------------
   # This solver for a system of ordindary differential equations
   # was written for didactical purposes. It does not represent 
   # the state of the art in numerical integration techniques
   # but should rather demonstrate how the simplest integration
   # technique can be implemented for relatively simple use.
   # Check the package "odesolve" available from http://www.r-project.org
   # for professional solvers for ordinary differential equations.
   #
   # Arguments:
   # rhs:       function returning the right hand side of the 
   #            system of differential equations as a function
   #            of the arguments x (state variables), t (time),
   #            and param (model parameters).
   # x.ini:     start values of the state variables.
   # param:     model parameters (transferred to rhs).
   # t.out:     set of points in time at which the solution 
   #            should be calculated; the solution at the first
   #            value in t is set to x.ini, the solution at subsequent
   #            values is calculated by numerical integration.
   # dt.max:    maximum time step; if the difference between points
   #            in time at which output has to be provided, t, is
   #            larger than dt.max, then this output interval is 
   #            divided into step of at most dt.max internally to
   #            improve the accuracy of the solution at the next 
   #            output point.
   # algorithm: right now, the only options are "euler" (explicit first order
   #            Euler algorithm) and "euler.2.order" (explicit second order
   #            Euler algorithm).
   #
   # Return Value:
   # matrix with results for state variables (columns) at all output
   # time steps t.out (rows).
   #
   #                                        Peter Reichert    Dec.  22, 2002
   #                                        last modification April 09, 2006
   # -----------------------------------------------------------------------
 
   # determine number of equations and number of time steps:
   num.eq <- length(x.ini)
   steps <- length(t.out)

   # define and initialize result matrix:
   x <- matrix(nrow=steps,ncol=num.eq)
   colnames(x) <- names(x.ini)
   rownames(x) <- t.out
   x[1,] <- x.ini

   # perform integration:
   if ( algorithm == "euler" )
   {
      for ( i in 2:steps )
      {
         if ( is.na(dt.max) || dt.max >= t.out[i]-t.out[i-1] )
         {
            # output interval <= dt.max (or dt.max not available):
            # perform a single step to the next output time:
            x[i,] <- x[i-1,] + (t.out[i]-t.out[i-1])*rhs(x[i-1,],t.out[i-1],param,...)
         }
         else
         {
            # output interval > dt.max:
            # perform multiple steps of maximum size dt.max:
            x[i,] <- x[i-1,]
            steps.internal <- ceiling((t.out[i]-t.out[i-1])/dt.max)
            dt <- (t.out[i]-t.out[i-1])/steps.internal
            for ( j in 1:steps.internal )
            {
               t.current <- t.out[i-1] + (j-1)*dt
               x[i,] <- x[i,] + dt*rhs(x[i,],t.current,param,...)
            }
         }
      }
   }
   else
   {
      if ( algorithm == "euler.2.order" )
      {
         for ( i in 2:steps )
         {
            if ( is.na(dt.max) || dt.max >= t.out[i]-t.out[i-1] )
            {
               # output interval <= dt.max (or dt.max not available):
               # perform a single step to the next output time:
               t.mid <- 0.5*(t.out[i-1]+t.out[i])
               x.mid <- x[i-1,] + 0.5*(t.out[i]-t.out[i-1])*rhs(x[i-1,],t.out[i-1],param,...)
               x[i,] <- x[i-1,] + (t.out[i]-t.out[i-1])*rhs(x.mid,t.mid,param,...)
            }
            else
            {
               # output interval > dt.max:
               # perform multiple steps of maximum size dt.max:
               x[i,] <- x[i-1,]
               steps.internal <- ceiling((t.out[i]-t.out[i-1])/dt.max)
               dt <- (t.out[i]-t.out[i-1])/steps.internal
               for ( j in 1:steps.internal )
               {
                  t.current <- t.out[i-1] + (j-1)*dt
                  t.mid <- t.current + 0.5*dt
                  x.mid <- x[i,] + 0.5*dt*rhs(x[i,],t.current,param,...)
                  x[i,] <- x[i,] + dt*rhs(x.mid,t.mid,param,...)
               }
            }
         }
      }
      else
      {
         stop(paste("ode: algorithm \"",algorithm,"\" not implemented",sep=""))
      }
   }

   # return result matrix: 
   return(x)
}


# =========================================================================== #