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

_______________________________________________

R-sig-ecology mailing list

[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology