[R-sig-dyn-mod] ReacTran package - unexpected results
Julien Lavenus
julien.lavenus at ips.unibe.ch
Mon Nov 2 15:20:59 CET 2015
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
--
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_50[1,2:dim(Fick_solved)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5))
par(new=T)
plot(Fick_solved_50[2,2:dim(Fick_solved)[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)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5),col="red")
Fick_solved_50[2,2:dim(Fick_solved)[2]]==Fick_solved_20000[2,2:dim(Fick_solved)[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,100),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