連立常微分方程式

  • \frac{d \mathbf{x}}{dt} = M \mathbf{x}なる常微分方程式
  • M = V S V^{-1}と分解したうえで
  • \mathbf{x(t=t)} = V exp(St) V^{-1} \mathbf{x(t=0)}と表される
  • ただしexp(St)固有値\lambda_iについて、exp(\lambda_i t)を対角成分とする行列である
  • こんなに単純で良ければ、\mathbf{x}が時系列データとして観察されれば、時間差分を時間微分に代用して考えて、線形回帰すればよい
  • やってみる
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]]