Skip to content

Differential Equations Using R?

3 messages · nfruge@macalester.edu, Simon Wood, Brian Ripley

nfruge@macalester.edu
# 
To whom it may concern,
        I am a student at Macaleste College, and next semester Macalester
is going to offer a course for CellBio that is mainly statistically based.
For the most part the students will be using R for analysis. The problem is
there will be some simple differential equations for the students to solve.
The committee that in charge of the classes corriculam would like only to
have to use one software package "R" for the entire corse. Thus the
students will not have to learn another platform in order to solve
differential equations. I have been nominated to find out if there exists
such a program for "R". Hence the question: Does there exist a program,
which has been written for R, where differential equations can be solved
numerically? If there is please contact me at your convience on where I may
find this program. If there is not such a program would it be simpler to
write such a program or to just learn how to use MatLab or another math
program? Your input would be greatly appreciated. Thank you all for your
time.

Regards,
Nick Fruge


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Simon Wood
#
On Mon, 10 Sep 2001 nfruge at macalester.edu wrote:

            

      
      
      
    

        
For simple stuff I use this adaptive timestepping rk2(3) integrator:

ode<-function (s0, t0 = 0, t1 = 1, c = 0, odt = 0.1, tol = array(1e-05, 
    length(s0)), sup = ode.setup()) 
{
    if (sup$on) {
        if (sup$phase) 
            sup$points <- length(sup$i)/2
        else sup$points <- length(sup$i)
        if (sup$points != round(sup$points)) 
            stop("Something wrong with your plot index vector")
        if (sup$phase == FALSE) 
            sup$xl <- c(t0, t1)
        plot(sup$xl, sup$yl, type = "n", xlab = "", ylab = "")
    }
    a2 <- 0.25
    a3 <- 27/40
    b21 <- 0.25
    b31 <- -189/800
    b32 <- 729/800
    b41 <- 214/891
    b42 <- 1/33
    b43 <- 650/891
    cc1 <- 533/2106
    cc3 <- 800/1053
    cc4 <- -1/78
    s <- s0
    t <- t0
    k <- floor((t1 - t0)/odt) + 1
    n.states <- length(s0)
    ot <- array(0, k)
    os <- array(0, c(k, n.states))
    ot[1] <- t
    os[1, ] <- s
    op.i <- 2
    next.op.t <- odt
    dt <- 0.01
    g <- grad(s, t, c)
    while (t < t1) {
        bct1 <- b21 * dt
        k1 <- g
        new.s <- s + k1 * bct1
        k2 <- grad(new.s, t + dt * a2, c)
        bct1 <- b31 * dt
        bct2 <- b32 * dt
        new.s <- s + k1 * bct1 + k2 * bct2
        k3 <- grad(new.s, t + dt * a3, c)
        bct1 <- b41 * dt
        bct2 <- b42 * dt
        bct3 <- b43 * dt
        new.s <- s + k1 * bct1 + k2 * bct2 + k3 * bct3
        k4 <- grad(new.s, t + dt, c)
        bct1 <- cc1 * dt
        bct2 <- cc3 * dt
        bct3 <- cc4 * dt
        error <- s + bct1 * k1 + bct2 * k3 + bct3 * k4 - new.s
        error <- abs(error)
        if (sum(error > tol) > 0) {
            dt <- dt * (min(tol/error)^0.3333) * 0.9
            if (dt <= 0) 
                stop("Timestep shrunk to zero")
        }
        else {
            g <- k4
            s <- new.s
            t <- t + dt
            if (min(error) > 0) 
                dt <- min(0.9 * dt * min(tol/error)^0.3333, 5 * 
                  dt)
            else dt <- 5 * dt
            if (t >= next.op.t) {
                ot[op.i] <- t
                os[op.i, ] <- s
                if (sup$on) {
                  if (sup$phase) {
                    if (sup$line) 
                      for (i in 1:sup$points) {
                        kx <- sup$i[2 * i - 1]
                        ky <- sup$i[2 * i]
                        x <- c(os[op.i - 1, kx], s[kx])
                        y <- c(os[op.i - 1, ky], s[ky])
                        lines(x, y)
                      }
                    else for (i in 1:sup$points) {
                      kx <- sup$i[2 * i - 1]
                      ky <- sup$i[2 * i]
                      points(os[op.i - 1, kx], os[op.i - 1, ky], 
                        pch = 19, col = "white")
                      points(s[kx], s[ky], pch = 19, col = "red")
                    }
                  }
                  else {
                    x <- c(ot[op.i - 1], t)
                    for (i in 1:sup$points) {
                      ky <- sup$i[i]
                      y <- c(os[op.i - 1, ky], s[ky])
                      lines(x, y)
                    }
                  }
                }
                op.i <- op.i + 1
                next.op.t <- next.op.t + odt
            }
            if ((t + dt) > next.op.t) 
                dt <- next.op.t - t
            if ((t + dt) > t1) 
                dt <- t1 - t
        }
    }
    ot[k] <- t
    os[k, ] <- s
    r <- list(t = ot, s = os)
}


ode.setup<-function (on = FALSE, phase = FALSE, i = 1, line = TRUE, xl =
c(0, 1), yl = c(0, 1)) 
{
    op <- list(on = on, phase = phase, i = i, line = line, xl = xl, 
        yl = yl)
}


Here's an example (planetary motion):

   grad<-function(s,t,c)
     { g<-s
       r3<-(s[1]^2+s[2]^2)^1.5  #radius
       g[1]<-s[3]    # x
       g[2]<-s[4]    # y 
       g[3]<- -s[1]/r3     # dx/dt
       g[4]<- -s[2]/r3     # dy/dt
       g
     }
     s0<-c(1,0,0,0.6)

ode(s0,0,20,0,0.05,1e-6,sup=ode.setup(on=T,phase=T,i=c(1,2),line=T,xl=c(-1,1),yl=c(-1,1)))

Documentation:

Numerical solution of ordinary differential equations.

Description:

     Use RK2(3) to solve the system of equations whose gradient vector
     is given by grad(). The equations are of the form:

     ds_1/dt = g_1(s,t)

     ds_2/dt = g_2(s,t)

     etc.

     where s is the vector of state variables.

Usage:


ode(s0,t0=0,t1=1,c=0,odt=0.1,tol=array(1e-5,length(s0)),sup=ode.setup())

Arguments:

      s0: Initial values of state variables.

      t0: The start time for model integration.

      t1: The stop time for model integration. 

       c: Array of any parameters to be passed to grad().

     odt: The time interval at which to output data.

     tol: An array of tolerance values for the state variables (can use
          a single number)

     sup: A list controlling the display of the solutions

Details:

     This function solves the system using an adaptive rk2(3) scheme.
     The user needs to define a function `grad(s,c,t)' that provides
     the  gradient values defining the model (see below).

Value:

     The returned value is a list with elements `t' - a one dimensional
     array  of output times and  `s' - a 2 dimensional array of states.

Simon

  ______________________________________________________________________

            

      
      
      
    

        
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Brian Ripley
    

    

      
      #
    

  


  

      
On Mon, 10 Sep 2001 nfruge at macalester.edu wrote:

        

            

      
      
      
    

        
Have you looked at the odesolve package on CRAN?

auk% cat odesolve/DESCRIPTION
Package: odesolve
Version: 0.5-5
Date: 2001/05/09
Title: Solvers for Ordinary Differential Equations
Author: R. Woodrow Setzer <setzer.woodrow at epa.gov> with fixups by Martin
Maechler <maechler at stat.math.ethz.ch>
Maintainer: R. Woodrow Setzer <setzer.woodrow at epa.gov>
Depends: R (>= 1.0.0)
Description: This package provides an interface for the ODE solver lsoda.
        ODEs are expressed as R functions.
License: GPL version 2

Looks like it meets your needs.