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")