Hi Lilith, your script does exactly what is expected from an ODE system, i.e. the new state depends on the *past* state and the derivatives. You cannot set a new state directly to a fixed value in a "pure" ODE system. Fortunately, package deSolve has a solution for your problem, called "root finding" in combination with "events" -- see example 3 on the ?events help page. Hope it helps Thomas Petzoldt PS: The list address is r-sig-dynamic-models at r-project.org i.e. without "-owner". AN email to "-owners" is only sent to the list moderators.
On 10/31/2012 6:31 PM, Lilith wrote:
Hi There! I have tried to find the solution to this problem for a long
time and hope you can help.
I have the following function containing some odes:
|myfunction<- function(t, state, parameters) {
with(as.list(c(state, parameters)),{
if (X>20) { # this is an internal threshold!
Y<- 35000
dY<- 0
}else{
dY<- b* (Y-Z)
}
dX<- a*X^6 + Y*Z
dZ<- -X*Y+ c*Y- Z
# return the rate of change
list(c(dX, dY, dZ),Y,dY)
})
}|
Here are some results:
|library(deSolve)
parameters<- c(a= -8/3, b= -10, c= 28)
state<- c(X= 1, Y= 1, Z= 1)
times<- seq(0, 10, by = 0.1)
out <- ode(y= state, times= times, func= myfunction, parms= parameters)
out
time X Y Z Y dY
1 0.0 1.000000 1.000000 1.000000 1.000000 0.00000
2 0.1 1.104670 2.132728 4.470145 2.132728 23.37417
3 0.2 1.783117 6.598806 14.086158 6.598806 74.87351
4 0.3 2.620428 20.325966 42.957134 20.325966 226.31169
5 0.4 3.775597 60.969424 126.920014 60.969424 659.50590
6 0.5 5.358823 176.094907 358.726482 176.094907 1826.31575
7 0.6 7.460841 482.506706 953.270570 482.506706 4707.63864
8 0.7 10.122371 1230.831764 2330.599161 1230.831764 10997.67398
9 0.8 13.279052 2859.284114 5113.458479 2859.284114 22541.74365
10 0.9 16.711405 5912.675147 9823.406760 5912.675147 39107.31613
11 1.0 24.452867 10590.600567 16288.435139 35000.000000 0.00000
12 1.1 25.988924 10590.600567 23476.343542 35000.000000 0.00000
13 1.2 26.572411 10590.600567 26821.703961 35000.000000 0.00000
14 1.3 26.844240 10590.600567 28510.668725 35000.000000 0.00000
15 1.4 26.980647 10590.600567 29391.032472 35000.000000 0.00000
...
|
States Y (columns 3 and 5) should be the same, instead when the
threshold is active they become different. Is there a way to "force"
Y(column 3) to change according to the set threshold?
Many thanks for your help!
Dr. Thomas Petzoldt Technische Universitaet Dresden Faculty of Environmental Sciences Institute of Hydrobiology 01062 Dresden, Germany E-Mail: thomas.petzoldt at tu-dresden.de http://tu-dresden.de/Members/thomas.petzoldt