[R] Differential Equations Using R?

Simon Wood snw at mcs.st-and.ac.uk
Tue Sep 11 07:27:37 CEST 2001


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

> 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
> 
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

  ______________________________________________________________________
> Simon Wood  snw at st-and.ac.uk  http://www.ruwpa.st-and.ac.uk/simon.html
> The Mathematical Institute, North Haugh, St. Andrews, Fife KY16 9SS UK
> Direct telephone: (0)1334 463799          Indirect fax: (0)1334 463748 



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list