1

R の deSolve の早期終了は、ルート関数を使用し、イベント関数を提供しないことで達成できることを知っています。これにより、ルートが見つかったときに統合が終了します。ただし、この手順を使用することにより、根を求める機能を備えたソルバーの適用に制限されます。

実際、ルートの正確な位置が問題にならない問題を扱っています。状態変数に突然の変化を導入する必要がありますが、これが発生する正確な瞬間は重要ではありません。したがって、条件が満たされたときに統合を停止し、突然の変化を導入して新しい開始状態ベクトルを再計算し、統合を再開することができます。これにより、deSolve パッケージで利用可能な多くのソルバーのいずれかを使用できるという柔軟性が得られます。

これを行うための推奨される方法はありますか?

編集

次の単純化された例を考えてみましょう。表現されたシステムは、1 の一定速度で 1 次元を移動するオブジェクトです。オブジェクトは位置 x=0 から開始し、次元の正の方向に移動します。オブジェクトが原点から 10 以上の距離に達した場合、位置は x=10 の点に関して参照されるように、座標の原点の変更を実行することを目的としています。これは、この時点の位置から 10 を引くことで単純化できます。

ルートを使用すると、これは次のように実現できます。

library(deSolve)
odeModel <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        x <- state
        dx <- 1
        list(
            dx
        )
    })
}

initialState <- 0
times <- 0:15

rootFunc <- function(t, state, parameters) {
    return(state-10)
}

eventFunc <- function(t, state, parameters) {
    return(state - 10)
}

integrationTest <- ode(y=initialState, times=times, func=odeModel, parms=NULL, 
                       events=list(func=eventFunc, root=TRUE), rootfun=rootFunc)

ただし、上記の理由から、ルートを使用せずにこれを実行しようとしています。条件が満たされているかどうかを追跡し、積分の結果を調べるいくつかの変数を出力に含めることによってそれを行うことができます (数値でなければなりません。そうしないと、deSolve は各評価時間の出力としてそれを受け入れません)。条件がいつ満たされたかを特定し、目的の変更を適用し、その時点から統合をやり直し、出力を結合します。前と同じ例を使用します。

library(deSolve)
distanceHigherThan10 <- function(x) {
    result <- if(x >= 10) 1 else 0
    return(result)
}

odeModel <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        x <- state
        dx <- 1
        higherThanThreshold <- distanceHigherThan10(x)
        list(
            dx, higherThanThreshold
        )
    })
}

initialState <- 0
times <- 0:15

integration <- ode(y=initialState, times=times, func=odeModel, parms=NULL)

changePoint <- which(diff(integration[,3]) != 0)[1] + 1
newIntegrationTimes <- times[changePoint:length(times)] - times[changePoint]
newStartingPosition <- integration[changePoint, 2] - 10
newIntegration <- ode(y=newStartingPosition, times=newIntegrationTimes, func=odeModel, parms=NULL)

newIntegration[, 1] <- newIntegration[, 1] + times[changePoint]
newIntegration[, 3] <- 1

combinedResults <- rbind(integration[1:(changePoint-1), ], newIntegration)

ただし、time=10以降の積分を2回行うため無駄な計算になります。条件が満たされた後に最初の統合を停止し、2 番目の統合部分に直接進む方が効率的です。これが、特定の条件が満たされたときに統合を停止し、その時点までの結果を返す方法があるかどうかを尋ねている理由です。

4

1 に答える 1