rk4                 package:odesolve                 R Documentation

_S_o_l_v_e _S_y_s_t_e_m _o_f _O_D_E (_o_r_d_i_n_a_r_y _d_i_f_f_e_r_e_n_t_i_a_l _e_q_u_a_t_i_o_n)_s _b_y
_c_l_a_s_s_i_c_a_l _R_u_n_g_e-_K_u_t_t_a _4_t_h _o_r_d_e_r _i_n_t_e_g_r_a_t_i_o_n.

_D_e_s_c_r_i_p_t_i_o_n:

     Solving initial value problems for  systems of first-order
     ordinary differential equations (ODEs) using the classical
     Runge-Kutta 4th order integration. The system of ODE's may be
     written as an R function (which may, of course, use '.C',
     '.Fortran', '.Call', etc., to call foreign code). A vector of
     parameters is passed to the ODEs, so the solver may be used as
     part of a modeling package for ODEs, or for parameter estimation
     using any appropriate modeling tool for non-linear models in R
     such as 'optim', 'nls', 'nlm' or 'nlme'.

_U_s_a_g_e:

       rk4(y, times, func, parms)

_A_r_g_u_m_e_n_t_s:

       y: the initial values for the ode system.  If 'y' has a name
          attribute, the names will be used to label the output matrix.

   times: times at which explicit estimates for 'y' are desired.  The
          first value in 'times' must be the initial time.

    func: a user-supplied function that computes the values of the
          derivatives in the ode system (the _model defininition_) at
          time t. The user-supplied function  'func' must be called as:
          'yprime = func(t, y, parms)'.  't' is the current time point
          in the integration, 'y' is the current estimate of the
          variables in the ode system, and 'parms' is a vector of
          parameters (which may have a names attribute, desirable in a
          large system).

          The return value of func should be a list, whose first
          element is a vector containing the derivatives of 'y' with
          respect to 'time', and whose second element is a vector
          (possibly with a 'names' attribute) of global values that are
          required at each point in 'times'. 

   parms: any parameters used in 'func' that should be modifiable
          without rewriting the function.

_D_e_t_a_i_l_s:

     The method is implemented primarily for didactic purposes. Please
     use lsoda for your real work!

_V_a_l_u_e:

     A matrix with up to as many rows as elements in 'times' and as
     many columns as elements in 'y' plus the number of "global" values
     returned in the second element of the return from 'func', plus and
     additional column for the time value.  There will be a row for
     each element in 'times' If 'y' has a names attribute, it will be
     used to label the columns of the output value.

_A_u_t_h_o_r(_s):

     Thomas Petzoldt petzoldt@rcs.urz.tu-dresden.de

_R_e_f_e_r_e_n_c_e_s:

     Press, W. H.,  Teukolsky, S. A., Vetterling, W. T. and Flannery,
     B. P. (1992) Numerical Recipes in C. Cambridge University Press.

_S_e_e _A_l_s_o:

     'lsoda'

_E_x_a_m_p_l_e_s:

     ## A simple resource limited Lotka-Volterra-Model
     lvmodel <- function(t, x, parms) {
       s <- x[1] # substrate
       p <- x[2] # producer
       k <- x[3] # consumer
       with(as.list(parms),{
         import <- approx(signal$times, signal$import, t)$y
         ds <- import - b*s*p + g*k
         dp <- c*s*p  - d*k*p
         dk <- e*p*k  - f*k
         res<-c(ds, dp, dk)
         list(res)
       })
     }

     ## vector of timesteps
     times  <- seq(0, 100, length=101)

     ## external signal with rectangle impulse
     signal <- as.data.frame(list(times = times,
                                 import = rep(0,length(times))))

     signal$import[signal$times >= 10 & signal$times <=11] <- 0.2

     ## Parameters for steady state conditions
     parms  <- c(a=0.0, b=0.0, c=0.1, d=0.1, e=0.1, f=0.1, g=0.0)

     ## Start values for steady state
     y<-xstart <- c(s=1, p=1, k=1)

     ## Classical RK4 with fixed time step
     out1  <- as.data.frame(rk4(xstart, times, lvmodel, parms))

     ## LSODA (default step size)
     out2 <- as.data.frame(lsoda(xstart, times, lvmodel, parms))

     ## LSODA: with fixed maximum time step
     out3 <- as.data.frame(lsoda(xstart, times, lvmodel, parms, hmax=1))

     par(mfrow=c(2,2))
     plot (out1$time, out1$s,  type="l",   ylim=c(0,3))
     lines(out2$time, out2$s, col="red",   lty="dotted")
     lines(out3$time, out3$s, col="green", lty="dotted")

     plot (out1$time, out1$p, type="l",    ylim=c(0,3))
     lines(out2$time, out2$p, col="red",   lty="dotted")
     lines(out3$time, out3$p, col="green", lty="dotted")

     plot (out1$time, out1$k, type="l",    ylim=c(0,3))
     lines(out2$time, out2$k, col="red",   lty="dotted")
     lines(out3$time, out3$k, col="green", lty="dotted")

     plot (out1$p, out1$k, type="l")
     lines(out2$p, out2$k, col="red",   lty="dotted")
     lines(out3$p, out3$k, col="green", lty="dotted")

