differences in ode.1D and tran.1D

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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

_______________________________________________
R-sig-ecology mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Reply | Threaded
Open this post in threaded view
|

Re: differences in ode.1D and tran.1D

Thomas Petzoldt
Hi Jeremy,

I saw this question today and hope it is not too late. Here a few thoughts:

1) Accuracy of computers is limited and rounding errors are unavoidable,
double precision numbers have "15–17 significant decimal digits
precision" (cf. Wikipedia), i.e. roughly a relative precision of 10^-16,
see:

https://en.wikipedia.org/wiki/Double-precision_floating-point_format

The default precision of the solvers is atol=1e-6, rtol=1e-6. You may
consider to change this, but it has of course limits.

2) For advective-diffusion equations, it is fundamental to fulfil the
CFL (Courant–Friedrichs–Lewy) condition:

https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition

It is not internally checked.

3) Package "ReacTran" has a function function advection.1D with some
higher order schemes. They have of course pros and cons, that are widely
discussed in the literature. The example at ?advection.1D gives you a a
few impressions.

4) More details may be discussed at the "R-SIG-Dynamic-Models" mailing list:

https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

Best regards,

Thomas


Am 15.02.2017 um 18:09 schrieb 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
>


--
Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

E-Mail: [hidden email]
http://tu-dresden.de/Members/thomas.petzoldt

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