r - 如何使用内部阈值解决 ODE？

Question

我有以下包含一些颂歌的函数:

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)
    })
}

以下是一些结果:

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
...

状态Y不同，有人可以解释我为什么吗？

我相信我没有正确设置我的阈值。
有没有办法做到这一点？

谢谢!

Answer 1

考虑解决 ODE 的最简单方法，即 Euler 方法:

 state = state+myfunction(t,state,parameters)*h
 f(t+h)=f(t) + f'(t) *h 

h是一个小的时间步长， myfunction是 f'(t) f(t) 的衍生物并且只计算导数，它无权访问实际 state也不是 Y .两者都在 ode 内部设置使用原则上类似于欧拉的方法:给定 f(t),f'(t),h 的数值它只是更新状态 f(t+h) .

所以阈值调整 dY但无法访问 state["Y"] .该过程只是操作一个局部变量，它被评估为 35000在 dX <- a*X^6 + Y*Z和 dZ <- -X*Y + c*Y - Z但实际 state["Y"]在 myfuction 之后被覆盖已返回 ode 内功能。

恐怕我想不出一个简单的方法来绕过这个设计。我只会用 out[5] .