2016-06-20 8 views
3

matlabでODE45などの可変時間ステップソルバーを使用している場合、出力の時間幅を定義することになります。つまり、times = [0 50]となります。ODE Times Matlab vs R

しかし、RではODEに結果を返す時間点を定義する必要があります。つまり、times = 0:50を指定すると、0,1,2, ... 50に51個の結果が返されます。他の賢明な私は、times = seq(0,50,0.1)のようなシーケンスを供給する必要があります。

私は、最初は急速に変化し、その後は徐々に変化する機能を持っています。 MATLABでは、出力結果には結果に82のタイムステップが返され、そのうちの49がタイムステップ0と1の間にあります。

結果を返す方法があるかどうかを知りたいMATLABと同じように、私は時間指定をあらかじめ指定しなくても、結果は返されたかったのです。

答えて

2

R - パッケージdeSolveは、このように使用されるパラメータtimesを説明:出力が望まれる

時系列います。

デニスは、もう一つの重要な文があるのが好きドキュメント:

我々はいくつかの実装のために、出力が望まれるれる ベクトルがメソッドれるメッシュを定義する、ということに注意してくださいに そのステップを実行するので、解の精度は入力ベクトルの時間に大きく依存します。

簡単な例は、パッケージdeSolveについての本で述べた以下の1、ローレンツ方程式、次のとおりです。

library(deSolve) 

parameters <- c(
    a = -8/3, 
    b = -10, 
    c = 28) 

state <- c(
    X = 1, 
    Y = 1, 
    Z = 1) 

# ---- define function in R 
Lorenz <- function(t, state, parameters) { 
    with(as.list(c(state, parameters)),{ 
    # rate of change 
    dX <- a*X + Y*Z 
    dY <- b * (Y-Z) 
    dZ <- -X*Y + c*Y - Z 

    # return the rate of change 
    list(c(dX, dY, dZ)) 
    }) # end with(as.list ... 
} 

times_1 <- seq(0, 100, by = 1) 
out_1 <- lsoda(y = state, times = times_1, func = Lorenz, parms = parameters) 

times_2 <- seq(0, 100, by = 0.01) 
out_2 <- lsoda(y = state, times = times_2, func = Lorenz, parms = parameters) 

tail(out_1) 

     time  X   Y   Z 
[96,] 95 30.54833 11.802554 12.445819 
[97,] 96 21.26423 4.341405 4.785116 
[98,] 97 33.05220 13.021730 12.934761 
[99,] 98 21.06394 2.290509 1.717839 
[100,] 99 10.34613 1.242556 2.238154 
[101,] 100 32.87323 -13.054632 -13.194377 

tail(out_2) 

      time  X   Y   Z 
[9996,] 99.95 17.16735 -7.509679 -12.08159 
[9997,] 99.96 17.66567 -7.978368 -12.77713 
[9998,] 99.97 18.26620 -8.468668 -13.47134 
[9999,] 99.98 18.97496 -8.977816 -14.15177 
[10000,] 99.99 19.79639 -9.501998 -14.80335 
[10001,] 100.00 20.73260 -10.036203 -15.40840 

あなたはこれが来るトン= 100で結果の違いを見ることができます異なる定義times

ありがとうございます。
J_F

1

DESOLVEにthis documentを読み込むと、それは述べて:

我々はすべての0.01日出力を生成し、100日間IVPを解きます。この 滑らかな図形を得るためには、小さな出力ステップが必要です。一般に これは統合のタイムステップに影響しません。これは通常、あなたが実際に入力としてtimes = seq(0,50,0.1)を使用する必要があるように思わそのようソルバー

によって決定 です。グラフ内の「面白い」点だけを表示したい場合は、実際のソルバー出力の出力を後処理する小さな関数を書くことができると思います。