[R-sig-eco] 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: thomas.petzoldt at tu-dresden.de
http://tu-dresden.de/Members/thomas.petzoldt