2016-09-25 3 views
1

時間tと空間xに渡る変数y(t、x)の進化を記述するpdeがあり、その進化を3次元図にプロットしたいとしますt、x、y)である。デソルブでは私はpdeを解決することができますが、私はこの種のダイアグラムをどのように取得するかについては考えていません。deSolve - 時間発展のプロットをプロットする

DESOLVEパッケージ命令の例は、、...、200であり、x = 1、yはアブラムシ、T = 0であり、以下、である...、60:

library(deSolve) 

Aphid <- function(t, APHIDS, parameters) { 
deltax <- c (0.5, rep(1, numboxes - 1), 0.5) 
Flux <- -D * diff(c(0, APHIDS, 0))/deltax 
dAPHIDS <- -diff(Flux)/delx + APHIDS * r 
list(dAPHIDS) 
} 

D <- 0.3 # m2/day diffusion rate 
r <- 0.01 # /day net growth rate 
delx <- 1 # m thickness of boxes 
numboxes <- 60 
Distance <- seq(from = 0.5, by = delx, length.out = numboxes) 

APHIDS <- rep(0, times = numboxes) 
APHIDS[30:31] <- 1 
state <- c(APHIDS = APHIDS) # initialise state variables 

times <-seq(0, 200, by = 1) 
out <- ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid") 

"アウト" t、y(x1)、y(x2)、... y(x60)のすべてのデータを含む行列を生成します。 (t、x)のyの進化と変動を示すサーフェスプロットを作成するにはどうすればよいですか?

答えて

0

パッケージを使用することによって方法が少し変わります。しかし、out[,-1]はサーフェイスを描画するのに理想的なマトリックスフォームなので、わずかなコストで行うことができます。私はrglplot3Dパッケージを使用して2つの例を示しました。

out2 <- out[,-1] 
AphID <- 1:ncol(out2) 

library(rgl) 
persp3d(times, AphID, out2, col="gray50", zlab="y") 
# If you want to change color with value of Z-axis 
# persp3d(times, AphID, out2, zlab="y", col=topo.colors(256)[cut(c(out2), 256)]) 

library(plot3D) 
mat <- mesh(times, AphID) 
surf3D(mat$x, mat$y, out2, bty="f", ticktype="detailed", xlab="times", ylab="AphID", zlab="y") 

enter image description here

関連する問題