[R-sig-eco] differences in ode.1D and tran.1D

[Jeremy Chacon]

Hello excellent list,

I am trying to figure out how to use ReacTran to do some reaction-diffusion
modeling.  I am starting by implementing the examples. I am finding that
modeling diffusion using the example in ode.1D gives *slightly *different
results than using tran.1D, and would like to know which is more accurate,
or if I can do something to improve numerical accuracy.

As a side (and less important) question, I have previously implemented some
simple finite-differences methods for solving diffusion equations, and have
always had to make sure step-sizes were appropriate for stability and
accuracy.  Does anyone know whether this is done internally, and how?

Here is some code for modeling exponential growth + diffusion with or
without using tran.1D, showing that the numbers are a bit different in each
case.

library(ReacTran)

# diffusion + exponential growth in 1D, using ode.1d (note that here, I'm
specifying the diffusion)
# taken verbatim from the deSolve docs
Aphid = function(t, APHIDS, parms){
  deltax = c(0.5, rep(1, numboxes-1), 0.5)
  Flux = -D * diff(c(0,APHIDS,0)) / deltax
  dAPHIDS = -diff(Flux) / delx + APHIDS * r
  list(dAPHIDS)
}

D = 0.3 # m2 / day (diffusion rate)
r = 0.01 # / day (growth rate)
delx = 1 # m (box size)
numboxes = 60
# the spatial grid, "distance of boxes on plant, m, 1 m intervals" (only
gets used for plotting I think)
Distance = seq(from = 0.5, by = delx, length.out = numboxes)

# init conditions
APHIDS = rep(0, times = numboxes)
APHIDS[30:31] = 1
state = c(APHIDS = APHIDS)

# run
times = seq(0, 200, by = 1)
out1 = ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid")
image(out1, method = "filled.contour", grid = Distance)

# now let's do the same  thing with reacTran

aphid_model = function(t, aphid, parms = NULL){
  diffusion = tran.1D(C = aphid, flux.up = NULL, flux.down = NULL,
                      D = 0.3, dx = 1)$dC
  growth = aphid * 0.01
  list(dCdt = diffusion + growth)
}



aphid = rep(0, 60)
aphid[30:31] = 1
state = c(aphid = aphid)
times = seq(0,200, by = 1)
out2 = ode.1D(state, times, aphid_model, parms = 0, nspec = 1)
image(out2, method = "filled.contour")
all((out1 - out2) < 1e-16) ## same results with this accuracy
all((out1 - out2) < 1e-17) ## diff results with this accuracy

-- 

*___________________________________________________________________________Jeremy
M. Chacon, Ph.D.*

*Post-Doctoral Associate, Harcombe Lab*
*University of Minnesota*
*Ecology, Evolution and Behavior*

	[[alternative HTML version deleted]]