[R-sig-dyn-mod] ReacTran package - unexpected results
Julien Lavenus
julien.lavenus at ips.unibe.ch
Mon Nov 2 17:13:04 CET 2015
Dear all,
My sincere apologies, I realised I sent you a version of the code with
mistakes. Please find attached a corrected R script.
Many thanks in advance for your help,
Best regards,
Julien Lavenus
On 02/11/15 15:20, Julien Lavenus wrote:
> Dear all,
>
> I am trying to use the ReacTran package to solve a system of partial
> differential equations. To make sure I was using the package
> correctly, I first implemented two simple 1D models (see below and
> code attached). I took as initial conditions for the molecule
> concentration a centered normal distribution for both models.
> Unexpectedly, none of the two test-models gave the expected results:
>
> The first model had diffusion (D=1000) and no advection (v=0). I
> expected the distribution of the modeled entity to get more and more
> uniform as the time increased, and in the end to get perfectly flat.
> However for a reason, I do not understand, the system reaches a steady
> states showing a peak (The distribution at t=20 and t=20000 are
> strictly equal).
>
> In the second test, I set the diffusion to 0 and I took a constant
> positive advection velocity, with zero-concentration boundary
> conditions. I was expecting the peak not to change shape and to move
> right until it gets absorbed by the boundary. Here again, I obtained a
> steady state. The peak moves right (and flattens a bit) and then stays
> where it is.
>
> If you have any idea of what I could be doing wrong, I would be really
> grateful for your help,
> Many thanks in advance,
>
> Best regards,
>
> Julien Lavenus
>
>
>
>
>
>
>
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
--
Dr. Julien Lavenus
Post-doc in Pr. Cris Kuhlemeier's lab
University of Bern
Institute of Plant Sciences
Alterbergrain 21
3006 Bern
Switzerland
Tel: + 41 31 631 4965
-------------- next part --------------
###
install.packages("ReacTran")
library(ReacTran)
#####################################################
## FIRST SIMPLEST MODEL: Fick second law for diffusion
# d^2(C(x,t))/dx^2 = D. d^2(C(x,t))/dx^2
# x: position
# t: time
# D: Diffusion coefficient
# C(x,t): Concentration at position x and time t
# Model definition
Fick_model=function(t=0,y,parms=NULL)
{
tran= tran.1D(C=y,D=1000,dx=grid)$dC
return(list(dC=tran))
}
# Grid definition
C=dnorm(seq(-10,10,.1))
L <- 200 # length of the diffusion space
N <- length(C) # number of grid layers
grid <- setup.grid.1D(x.up = 0,L = L, N = N)
# Initial conditions + simulation
Fick_solved_20 <- ode.1D(y = C, times=seq(0,20),func=Fick_model,nspec=1,method="lsoda")
Fick_solved_20000 <- ode.1D(y = C, times=seq(0,20000),func=Fick_model,nspec=1,method="lsoda")
# Plot Results
plot(Fick_solved_20[1,2:dim(Fick_solved_20)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5))
par(new=T)
plot(Fick_solved_20[2,2:dim(Fick_solved_20)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5),col="blue")
par(new=T)
plot(Fick_solved_20000[2,2:dim(Fick_solved_20000)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5),col="red")
Fick_solved_20[2,2:dim(Fick_solved_20)[2]]==Fick_solved_20000[2,2:dim(Fick_solved_20000)[2]]
###
# Model definition
Advection_model=function(t=0,y,parms=NULL)
{
tran= tran.1D(C=y,D=0,dx=grid,v=10)$dC
return(list(dC=tran))
}
# Grid definition
C=dnorm(seq(-10,10,.1))
L <- 201 # length of the diffusion space
N <- length(C) # number of grid layers
grid <- setup.grid.1D(x.up = 0,L = L, N = N)
# Initial conditions + simulation
Advection_model_solved <- ode.1D(y = C, times=seq(0,1000),func=Advection_model,nspec=1,method="lsoda",)
# Plot Results
plot(Advection_model_solved[1,2:dim(Advection_model_solved)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5))
par(new=T)
plot(Advection_model_solved[2,2:dim(Advection_model_solved)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5),col="blue")
More information about the R-sig-dynamic-models
mailing list