連立常微分方程式

• なる常微分方程式は
• と分解したうえで
• と表される
• ただしは固有値について、を対角成分とする行列である
• こんなに単純で良ければ、が時系列データとして観察されれば、時間差分を時間微分に代用して考えて、線形回帰すればよい
• やってみる

d <- 4
A <- matrix(rnorm(d^2),d,d)
eigen.out <- eigen(A)
V <- eigen.out[[2]]
V.inv <- solve(V)
x.init <- rnorm(d)
t <- seq(from=-3,to=3,length=100)
X <- matrix(0,length(t),d)
for(i in 1:length(t)){
 X[i,] <- Re(V %*% diag(exp(eigen.out[[1]]*t[i])) %*% V.inv %*% x.init)
}
dX <- apply(X,2,diff)
lm.out <- lm(dX~X[-1,]-1) # 切片項は入れない
# もしくは
X.mean <- (X[-length(t),]+X[-1,])/2 # こちらの方が推定はよりよい
lm.out.2 <- lm(dX~X.mean-1)
lm.out
lm.out.2
eigen.out.post <- eigen(lm.out$coefficients)
eigen.out.post.2 <- eigen(lm.out.2$coefficients)
eigen.out.post
eigen.out.post.2
eigen.out
eigen.out[[1]] * (t[2]-t[1])
eigen.out.post[[1]]
eigen.out.post.2[[1]]