3

可変ステップ サイズの ODE 問題を解くために、 deSolve R パッケージの明示的なルンゲクッタ法ode45 (別名 rk45dp7)を使用したいと考えています。

deSolve のドキュメントによると、等距離時間ステップの代わりに ode45 メソッドを使用してrkソルバー関数に適応または可変時間ステップを使用することは可能ですが、これを行う方法がわかりません。

rk 関数は次のように呼び出されます。

rk(y, times, func, parms, rtol = 1e-6, atol = 1e-6, verbose = FALSE, tcrit = NULL,
hmin = 0, hmax = NULL, hini = hmax, ynames = TRUE, method = rkMethod("rk45dp7", ... ), 
maxsteps = 5000, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, 
ipar = NULL, nout = 0, outnames = NULL, forcings = NULL, initforc = NULL, fcontrol = 
NULL, events = NULL, ...)

timesは、y の明示的な推定が必要な時間のベクトルです。

距離が 0.01 の等間隔の時間ステップの場合、時間は次のように記述できます。

times <- seq(0, 100, 0.01)

0 から 100 までの間隔の方程式を解きたいとすると、ステップ サイズを指定せずに時間を定義するにはどうすればよいでしょうか?

どんな助けでも大歓迎です。

4

1 に答える 1

5

ここには 2 つの問題があります。最初に、複数のインクリメントを持つ時間のベクトルを指定する場合は、これを使用します (例):

times <- c(seq(0,0.9,0.1),seq(1,10,1))
times
#  [1]  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0 10.0

ここで、[0,1] は 0.1 で、[1,10] は 1 です。

ただし、実際にこれを行う必要はありません。パラメーターtimes=はrk(...)、結果を報告する時間を示します。適応アルゴリズムは、時間の増分を内部的に調整して、引数で指定された時間に正確な結果を生成します。したがって、適応アルゴリズムの場合、たとえば、method="rk45dp7"何もする必要はありません。たとえば、非適応アルゴリズムの場合、method="euler"アルゴリズムで使用される時間増分は、実際には で指定された増分times=です。この効果は、Van der Pol 発振器を統合したこの単純な例で確認できます。

y.prime <- function(t,y.vector,b) {    # Van der Pol oscillator
  x <- y.vector[1]
  y <- y.vector[2]
  x.prime <- y
  y.prime <- b*y*(1-x)^2 - x
  return(list(c(x=x.prime,y=y.prime)))
}
h  <- .001                  # time increment
t  <-  seq(0,10,h)          # times to report results
y0 <- c(0.01,0.01)          # initial conditions
euler   <- rk(y0, t,func=y.prime,parms=1,method="euler")
rk45dp7 <- rk(y0, t,func=y.prime,parms=1, method="rk45dp7")
# plot x vs. y
par(mfrow=c(1,2))
plot(euler[,2],euler[,3], type="l",xlab="X",ylab="Y",main=paste("Euler: h =",format(h,digits=3)))
plot(rk45dp7[,2],rk45dp7[,3], type="l",xlab="X",ylab="Y",main=paste("RK45DP7: h =",format(h,digits=3)))

以下は、 のいくつかの値に対する結果の比較ですh。method="rk45dp7"で結果が安定していることに注意してくださいh < 0.5。これは、rk45dp7が必要に応じて時間の増分を内部で調整しているためです。結果はまでmethod="euler"一致しません。rk45dp7h~0.01

于 2013-12-15T16:12:47.910 に答える