混乱するかもしれませんが、次のように完全にうまく処理できる連立方程式を説明しているようです( ODElsoda
を実装しましたが、何を考えているのかわからなかったため、いくつかのパラメーターを作成しました)。
gfun <- function(t,y,parms,...) {
## 'with' trick lets us write gradient in terms of variable/parameter names
with(as.list(c(y,parms)),
list(c(b=beta-k*b,a=alpha-b*gamma),NULL))
}
library(deSolve)
L1 <- lsoda(y=c(b=1,a=1),
times=seq(0,10,by=0.1),
func=gfun,
parms=c(alpha=0.1,beta=0.2,gamma=0.05,k=0.01))
matplot(L1[,1],L1[,-1],type="l",lty=1,bty="l",las=1)
PS:これは結合された線形常微分方程式のセットのようです。したがって、数値的に解くのではなく、実際には完全な閉形式の解を得ることができるはずです。(私は今それをするのが面倒です; b(t)はすぐに解くことができます(「アフィン」方程式)、a(t)は積分によって解くことができます。)