[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