[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