inst/tests/runit.ThreepSerialFeedback_linear_vs_nonlinear.R In SoilR: Models of Soil Organic Matter Decomposition

```#
# vim:set ff=unix expandtab ts=2 sw=2:
# This test usese the non linear model approach for a linear problem to show that the results are consistent
test.ThreepSerialFeedback_linear_vs_nonlinear=function(){
require(RUnit)
t_start=0
t_end=20
tn=100
tol=.02/tn
print(tol)
timestep=(t_end-t_start)/tn
t=seq(t_start,t_end,timestep)
k1=1/2
k2=1/3
k3=1
a12=1/9
a23=1/6
a21=1/9
a32=1/6
A=new("ConstLinDecompOp",matrix(
byrow=T,
nrow=3,
ncol=3,
c(
-k1,   a12,    0,
a21,   -k2,  a23,
0,    a32,  -k3
)
)
)

alpha=list()
alpha[["2_to_1"]]=function(C,t){
a12/k2
}
alpha[["3_to_2"]]=function(C,t){
a23/k3
}
alpha[["1_to_2"]]=function(C,t){
a21/k1
}
alpha[["2_to_3"]]=function(C,t){
a32/k2
}
nr=3
N=matrix(
nrow=nr,
ncol=nr,
c(
k1,    0,     0,
0  ,  k2,     0,
0,     0,    k3
)
)
f=function(C,t){
# in this case the application of f can be expressed by a matrix multiplication
# f(C,t)=N C
# furthermorde the matrix N is actually completely linear and even constant
# so we can write f(C,t)  as a Matrix product
# note however that we could anything we like with the components
# of C here.
# The only thing to take care of is that we release a vector of the same
# size as C
return(N%*%C)
}
Anl=new("TransportDecompositionOperator",t_start,Inf,nr,alpha,f)

c01=3
c02=2
c03=1
iv=c(c01,c02,c03)
inputrates=new("TimeMap",t_start,t_end,function(t){return(matrix(
nrow=nr,
ncol=1,
c(
0,  2,  2
)
))})
#################################################################################
# we check if we can reproduce the linear decomposition operator from the
# nonlinear one
Tr=getTransferMatrix(Anl) #this is a function of C and t
T_00=Tr(matrix(nrow=nr,iv),0)
af=getFunctionDefinition(A)
af_0=af(0)
#pp("T_00",environment())
#pp("af_0",environment())
#pe(quote(T_00%*%N),environment())
checkEquals(af_0,T_00%*%N)

#################################################################################
# build the two models (linear and nonlinear)
mod=GeneralModel( t, A,iv, inputrates, deSolve.lsoda.wrapper)
modnl=GeneralNlModel( t, Anl, iv, inputrates, deSolve.lsoda.wrapper)
# compare the Cstock
Y=getC(mod)
Ynonlin=getC(modnl)
#R=getReleaseFlux(mod)
#Rnonlin=getReleaseFlux(modnl)
#begin plots
lt1=2
lt2=4
m=matrix(c(1,1),1,1,byrow=TRUE)
ex=expression(
layout(m),
plot(t,Y[,1],type="l",lty=lt1,col=1,ylab="Concentrations",xlab="Time",ylim=c(min(Y),max(Y))),
lines(t,Ynonlin[,1],type="l",lty=lt2,col=1),
lines(t,Y[,2],type="l",lty=lt1,col=2),
lines(t,Ynonlin[,2],type="l",lty=lt2,col=2),
lines(t,Y[,3],type="l",lty=lt1,col=3),
lines(t,Ynonlin[,3],type="l",lty=lt2,col=3),
legend(
"topright",
c(
"linear sol for pool 1",
"non linear sol for pool 1",
"linear sol for pool 2",
"non linear sol for pool 2",
"linear sol for pool 3",
"non linear sol for pool 3"
),
lty=c(lt1,lt2),
col=c(1,1,2,2,3,3)
)
)
plotAndCheck("runit.ThreepSerialFeedback_linear_vs_nonlinear.pdf",ex,environment())
#plot(t,R[,1],type="l",lty=lt1,col=1,ylab="Respirationfluxes",xlab="Time",ylim=c(min(R),max(R)))
#lines(t,Rnonlin[,1],type="l",lty=lt2,col=1)
#lines(t,R[,2],type="l",lty=lt1,col=2)
#lines(t,Rnonlin[,2],type="l",lty=lt2,col=2)
#lines(t,R[,3],type="l",lty=lt1,col=3)
#lines(t,Rnonlin[,3],type="l",lty=lt2,col=3)
#legend(
#"topright",
#  c(
#  "linear sol for pool 1",
#  "non linear sol for pool 1",
#  "linear sol for pool 2",
#  "non linear sol for pool 2",
#  "linear sol for pool 3",
#  "non linear sol for pool 3"
#  ),
#  lty=c(lt1,lt2),
#  col=c(1,1,2,2,3,3)
#)
# end plots
# begin checks
checkEquals(
Y,
Ynonlin,
"test non linear solution for C-Content computed by the ode mehtod against analytical",
tolerance = tol,
)
#checkEquals(
# R,
# Rnonlin,
# "test non linear solution for Respiration computed by the ode mehtod against analytical",
# tolerance = tol,
#)

}
```

Try the SoilR package in your browser

Any scripts or data that you put into this service are public.

SoilR documentation built on May 29, 2017, 10:57 a.m.