如何使用內部閾值來解決ODE？

我已經包含了一些賦有以下功能：如何使用內部閾值來解決ODE？

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是不同的，任何人都可以解釋我爲什麼嗎？

我相信我沒有正確設定我的門檻。有沒有辦法呢？

謝謝！

來源

2012-10-31 Claudia

我試圖破譯'ode'的幫助文件，但沒有成功。我所能建議的只是嘗試一個非常簡單的函數，並查看各種返回的列表示的內容。我懷疑這兩個「Y」列正在看不同的東西。 –

我修改了代碼以突出顯示第3列和第5列應該看起來相同。但是，當閾值變爲活動狀態（從第11行開始）時，它們會返回不同的結果？！ – Claudia

再一次，我不知道底層的'soda *'算法是幹什麼的，但我現在的猜測是'10590.600567'是前一個循環的輸出（當'dY'還不是零時），並且該值保留在任何「Y輸入」欄表示，而「Y輸出」正確地凍結在35000。 –

想到的最簡單的方法來解決常微分方程，即歐拉方法：

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"]。該過程僅處理一個局部變量，該局部變量在dX <- a*X^6 + Y*Z和dZ <- -X*Y + c*Y - Z中評估爲35000，但在myfuction已返回到ode函數內後，實際state["Y"]被覆蓋。

恐怕我想不出一種簡單的方法來繞過這個設計。我只會使用out[5]。

來源

2013-07-08 19:49:52 user2561915

如何使用內部閾值來解決ODE？

回答

相關問題