# [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
> 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 <- t
os[1, ] <- s
op.i <- 2
next.op.t <- odt
dt <- 0.01
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):

{ g<-s
g<-s    # x
g<-s    # y
g<- -s/r3     # dx/dt
g<- -s/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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

```