[R-sig-dyn-mod] a simulation using base R code
YA
x|nx|813 @end|ng |rom 126@com
Tue Jul 14 11:42:22 CEST 2020
hi Tom,
Thank you very much for the information. The output using deSolve is indeed a lot nicer, however, I still have not get what I expected. So as the hehavior increases it consumes more resource, as the resource decrease the increase of behavior should be slow down, which results in a flat curve of the behavior. Namely, the resource plays as a limiting factor for the growth of the behavior, so that the behavior would not be increasing exponentially forever.
My new code looks like below:
# model 1
library(deSolve)
simtime=seq(1,1000,by=1/8)
stocks=c(behavior=1,resource=5)
auxs=c(inprob.beh=0.01,outprob.beh=0,inprob.res=0,outprob.res=0.05)
myfun=function(time,stocks,auxs){
with(as.list(c(stocks,auxs)),{
resource=sqrt(behavior);
inf.beh=behavior*inprob.beh; outf.beh=behavior*outprob.beh; netf.beh=inf.beh-outf.beh;
inf.res=resource*inprob.res; outf.res=resource*outprob.res; netf.res=inf.res-outf.res;
list(c(netf.beh,netf.res),inflow.behavior=inf.beh,outflow.resource=outf.res)
})
}
out=ode(y=stocks,parms = auxs,time=simtime,func=myfun)
out
plot(out)
To connect the relationship between the behavior and resource, I used 'resource=sqrt(behavior);' I saw reference that sqrt() can be helpful for limiting growth, but it seems not working in my situation. Any other moves?
A successful code looks like this:
# model3
library(deSolve)
simtime=seq(0,100,by=1/8)
stocks=c(s.machine=100,s.econ=50)
auxs=c(a.labor=100)
myfun=function(time,stocks,auxs){
with(as.list(c(stocks,auxs)),{
s.econ=a.labor*sqrt(s.machine); # connect two stock variable into a relationship
inflow.machine=s.econ*0.2;outflow.machine=s.machine*0.1;netflow.machine=inflow.machine-outflow.machine;
inflow.econ=0.12;outflow.econ=0.04;netflow.econ=inflow.econ-outflow.econ;
list(c(netflow.machine,netflow.econ),inflow.machine,outflow.machine)
})
}
out=ode(y=stocks,parms = auxs,times=simtime,func=myfun)
out
plot(out)
Thank you very much.
best regards,
YA
------------------ Original ------------------
From: "Special Interest Group for Dynamic Simulation Models in R" <tom.shatwell using ufz.de>;
Date: Sun, Jul 12, 2020 08:29 PM
To: "r-sig-dynamic-models"<r-sig-dynamic-models using r-project.org>;
Subject: Re: [R-sig-dyn-mod] a simulation using base R code
Hi JA,
the problem might be that you are adding the two terms in the inflow.
There should be some multiplicative term to make them dependent. Try
simulating longer than 10 time units too.
You might like to try using the deSolve package. It's much better than
just using differencing loops. I didn't really get exactly what you are
trying to model, but there is some code below that would do something
similar - just adapt as necessary.
Cheers, Tom
# parameters
pars <- c(p1=0.02, p2=0.01, p3=-0.05)
# model
model <- function(t, state, parms) {
with(as.list(c(state, parms)), {
# rates
flow.stock <- stock*p1+resource*p2 # inflow
flow.resource <- stock*p3 # outflow
# derivs
dstock <- flow.stock * stock
dresource <- flow.resource*resource
list(c(dstock, dresource), c(flow.stock=flow.stock,
flow.resource=flow.resource))
})
}
# initial state
y <- c(stock=1, resource=5)
# times
unit <- 1/8
times <- seq(0,10,by=unit)
# simulation
library(deSolve)
out <- ode(y = y, times = times, func = model, parms = pars)
plot(out)
On 2020-07-12 12:07, YA wrote:
> hi Jingsong,
>
>
> Thanks for the quick response.
>
>
> How come the two loops are independent? I specified the following:
>
>
> > flow.stock=stock*0.02+resource*0.01 # inflow
> > flow.resource=stock*(-0.05) # outflow
>
>
> So that the flow.stock and flow.resource be influenced by each other. Why they didnt work?
>
>
> YA
>
>
>
>
>
>
>
> ------------------ Original ------------------
> From: "Jinsong Zhao" <jszhao using yeah.net>;
> Date: Sun, Jul 12, 2020 04:20 PM
> To: "YA"<xinxi813 using 126.com>;
>
> Subject: Re: [R-sig-dyn-mod] a simulation using base R code
>
>
>
> On 2020/7/12 15:56, YA wrote:
> > Dear list,
> >
> >
> > I am trying to simulate a system dynamics model with two state variables, stock and resource. I expect that as the resource draining out, the inflow of stock would decrease, then the increase of the stock would be slow down, but still keep increasing the outflow of the resource. So the system would eventually runing out of resource and the curve of the stock would getting flat. However, the simulation result is not consistent with my hypothesis. Could you give me some advice please? I used the base R code without any system dynamics package, is that ok? or do I need to use the function in the package to get the simulation results?
> >
> >
> > Here is a reproducible code:
> >
> >
> > stock=1
> > resource=5
> > flow.stock=stock*0.02+resource*0.01 # inflow
> > flow.resource=stock*(-0.05) # outflow
> > unit=1/8
> > resource=resource+flow.resource*resource
> > stock=stock+flow.stock*stock
> > for (i in 1:length(seq(1,10,by=unit))) stock[i+1]=stock[i]+flow.stock*stock[i]
>
> In the above line, flow.stock is not changed.
>
> > for (i in 1:length(seq(1,10,by=unit))) resource[i+1]=resource[i]+flow.resource*resource[i]
>
> In the above line, flow.resource is not changed.
>
> you use two independent loop to calculate each variable. however, in the
> system you hope to simulate, both are coupled.
>
> Please read some introduction materials about modeling in R, for
> example, Soetaert, K., and Herman, P.M.J. (2009). A practical guide to
> ecological modelling: using R as a simulation platform (Dordrecht:
> Springer).
>
> Best,
> Jinsong
>
> > stock=c(3,stock)
> > resource=c(5,resource)
> > dat=data.frame(id=1:length(stock),time=c(0,seq(1,10,by=unit),10+unit),stock,resource)
> > dat
> > plot(dat$stock)
> > lines(dat$stock)
> > lines(dat$resource)
> >
> >
> > Thank you very much.
> >
> >
> > Best regards,
> >
> >
> > YA
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-dynamic-models mailing list
> > R-sig-dynamic-models using r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> >
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
--
------------------------------------------------------------------------------------------
Dr. Tom Shatwell
Researcher
Department of Lake Research
Helmholtz-Zentrum für Umweltforschung GmbH - UFZ
Helmholtz Centre for Environmental Research GmbH - UFZ
Brueckstr. 3a, 39114, Magdeburg, Germany
Phone +49 391 810 9440, Fax +49 391 810 9150
tom.shatwell using ufz.de, www.ufz.de
Sitz der Gesellschaft/Registered Office: Leipzig
Registergericht/Registration Office: Amtsgericht Leipzig
Handelsregister Nr./Trade Register Nr.: B 4703
Vorsitzende des Aufsichtsrats/Chairwoman of the Supervisory Board: MinDirig'in Oda Keppler
Wissenschaftlicher Geschäftsführer/Scientific Managing Director: Prof. Dr. Georg Teutsch
Administrative Geschäftsführerin/Administrative Managing Director: Dr. Sabine König
Vermeiden Sie unnötige Ausdrucke./Think before printing.
_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models using r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
[[alternative HTML version deleted]]
More information about the R-sig-dynamic-models
mailing list