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 |
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 |
Free forum by Nabble | Edit this page |