[R-sig-dyn-mod] How to solve ODEs with an internal threshold?

[Thomas Petzoldt]

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