我有以下包含一些颂歌的功能:
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州不同，有人可以向我解释为什么吗？
我相信我没有正确设置阈值。
有办法吗？
谢谢!
最佳答案
考虑解决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中使用一种原则上类似于Euler的方法在内部设置:给定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]。关于r - 如何用内部阈值求解ODE？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/13154681/