---
title: "行列はありがたい"
author: "ryamada"
date: "2017年6月12日"
output:
html_document:
toc: true
toc_depth: 6
number_section: true
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
$$
\frac{\Delta x}{\Delta t} = p x + q y\\
\frac{\Delta y}{\Delta t} = r x + s y
$$
```{r}
delta.t <- 10^-3
n.step <- 10^4
p <- 1
q <- -2
r <- 3
s <- -2
M <- matrix(c(p,q,r,s),byrow=TRUE,ncol=2,nrow=2)
M
xy <- matrix(0,2,n.step)
xy[,1] <- c(3,5)
delta.xy <- M %*% xy[,1]
delta.xy
p * xy[1,1] + q * xy[2,1]
r * xy[1,1] + s * xy[2,1]
delta.xy * delta.t
```
```{r}
for(i in 2:n.step){
dxdy <- M %*% xy[,i-1] * delta.t
xy[,i] <- xy[,i-1] + dxdy
}
plot(xy[1,],xy[2,],type="l",asp=TRUE)
```
### 簡単に書ける、一般化できる
$$
\frac{\Delta x}{\Delta t} = p x + q y\\
\frac{\Delta y}{\Delta t} = r x + s y
$$
$$
\frac{\Delta x_1}{\Delta t} = m_{1,1} x_1 + m_{1,2} x_2 + m_{1,3} x_3\\
\frac{\Delta x_2}{\Delta t} = m_{2,1} x_1 + m_{2,2} x_2 + m_{2,3} x_3\\
\frac{\Delta x_3}{\Delta t} = m_{3,1} x_1 + m_{3,2} x_2 + m_{3,3} x_3
$$
$$
\frac{\Delta x_i}{\Delta t} = m_{i,1} x_1 + m_{i,2} x_2 + m_{i,3} x_3
$$
$$
\frac{\Delta x_i}{\Delta t} = \sum_{j=1}^3 m_{i,j} x_j
$$
$$
\begin{pmatrix}
\frac{\Delta x_1}{\Delta t}\\
...\\
\frac{\Delta x_n}{\Delta t}
\end{pmatrix} = \begin{pmatrix}
m_{1,1},m_{1,2},...,m_{i,n}\\
...\\
m_{n,1},m_{n,2},...,m_{n,n}
\end{pmatrix} \begin{pmatrix}
x_1\\
...\\
x_n
\end{pmatrix}
$$
$$
\frac{\Delta \mathbf{x}}{\Delta t} = M \mathbf{x}
$$
$x_1$と$x_2$の関係は2要素シミュレーションで用いた関係であり、
$x_3$と$x_4$の関係も同じだとする。
$x_2$と$x_3$とに「$x_2$は$x_3$の量に応じて増え」「$x_3$は$x_2$の量に応じて減る」という関係を入れてみよう。
2要素の行列
```{r}
delta.t <- 10^-3
n.step <- 10^4
p <- 1
q <- -2
r <- 3
s <- -2
M <- matrix(c(p,q,r,s),byrow=TRUE,ncol=2,nrow=2)
M
```
4要素に拡張
```{r}
M4 <- matrix(0,4,4)
M4[c(1,2),c(1,2)] <- M
M4
M4[c(3,4),c(3,4)] <- M
M4
# x2,x3の関係
a <- 3
b <- -2
M4[2,3] <- a
M4[3,2] <- b
M4
```
```{r}
x4 <- matrix(0,4,n.step)
x4[,1] <- c(3,5,3,5)
delta.x4 <- M4 %*% x4[,1]
delta.x4 * delta.t
```
```{r}
for(i in 2:n.step){
dx <- M4 %*% x4[,i-1] * delta.t
x4[,i] <- x4[,i-1] + dx
}
plot(x4[1,],x4[2,],type="l",asp=TRUE)
plot(x4[3,],x4[4,],type="l",asp=TRUE)
plot(x4[2,],x4[3,],type="l",asp=TRUE)
matplot(t(x4),type="l")
```
指数関数的に増える話。人口トピックのように。
$$
\frac{\Delta x}{\Delta t} = a x\\
a = \frac{x}{\frac{\Delta x}{\Delta t}}
$$
```{r}
a <- 0.2
n.step <- 10
x <- rep(0,n.step)
x[1] <- 1
for(i in 2:n.step){
delta.x <- a * x[i-1]
x[i] <- x[i-1] + delta.x
}
plot(x,type="b")
```
増える勢いは、隣合う値の割り算で求められた。
```{r}
x[-1]/x[-n.step]
a+1
```
微分(差分)方程式に照らしてこんな風にも出来る。
```{r}
diff(x) / x[-n.step]
a
```
$$
\frac{\Delta x}{\Delta t} = a x\\
a = \frac{\frac{\Delta x}{\Delta t}}{x}
$$
$$
\frac{\Delta \mathbf{x}}{\Delta t} = M \mathbf{x}\\
\begin{pmatrix}\frac{\Delta x}{\Delta t}\\ \frac{\Delta y}{\Delta t} \end{pmatrix} = M \begin{pmatrix} x \\ y \end{pmatrix}\\
M = \frac{\frac{\Delta \mathbf{x}}{\Delta t}}{\mathbf{x}}
$$
$$
\begin{pmatrix}\frac{\Delta x_1}{\Delta t}, \frac{\Delta x_2}{\Delta t}\\ \frac{\Delta y_1}{\Delta t}, \frac{\Delta y_2}{\Delta t} \end{pmatrix} = M \begin{pmatrix} x_1, x_2 \\ y_1, y_2 \end{pmatrix}\\
M = \frac{\begin{pmatrix}\frac{\Delta x_1}{\Delta t}, \frac{\Delta x_2}{\Delta t}\\ \frac{\Delta y_1}{\Delta t}, \frac{\Delta y_2}{\Delta t}\end{pmatrix}}{\begin{pmatrix} x_1, x_2 \\ y_1, y_2 \end{pmatrix}}\\
M = \begin{pmatrix}\frac{\Delta x_1}{\Delta t}, \frac{\Delta x_2}{\Delta t}\\ \frac{\Delta y_1}{\Delta t}, \frac{\Delta y_2}{\Delta t}\end{pmatrix} \times (\begin{pmatrix} x_1, x_2 \\ y_1, y_2 \end{pmatrix})^{-1}
$$
逆行列。
$$
\frac{1}{X} = X^{-1}
$$
```{r}
xy[,1:3]
xy[,1:2]
delta.xy12 <- xy[,2:3] - xy[,1:2]
xy12 <- xy[,1:2]
xy12.inv <- solve(xy12)
delta.xy12 %*% xy12.inv
M
M * delta.t
```
```{r}
M
```
```{r}
eigen.out <- eigen(M)
eigen.out
diag(eigen.out[[1]])
tmp <- eigen.out[[2]] %*% diag(eigen.out[[1]]) %*% solve(eigen.out[[2]])
round(tmp,5)
M
```
```{r}
M4
eigen.out <- eigen(M4)
eigen.out
diag(eigen.out[[1]])
tmp <- eigen.out[[2]] %*% diag(eigen.out[[1]]) %*% solve(eigen.out[[2]])
round(tmp,5)
M4
```
$$
M = V S U\\
U = V^{-1}\\
UV= VU = I\\
M^2 = M M = VSU VSU = VS(UV)SU = VSISU = VS^2U\\
M^k = VS^k U
$$
対角行列のk乗
```{r}
s <- eigen.out[[1]]
s
S <- diag(s)
S
```
```{r}
s2 <- s^2
s2
S2 <- diag(s^2)
S2
S %*% S
```
```{r}
k <- 5
Mout <- M4
for(i in 2:k){
Mout <- M4 %*% Mout
}
Mout
eigen.out[[2]] %*% diag(s^k) %*% solve(eigen.out[[2]])
```
$$
\frac{\Delta x}{\Delta t} = p x + q xy\\
\frac{\Delta y}{\Delta t} = r xy + s y
$$
```{r}
delta.t <- 10^-3
n.step <- 10^4
p <- 1
q <- -2
r <- 3
s <- -2
M <- matrix(c(p,q,r,s),byrow=TRUE,ncol=2,nrow=2)
M
xy <- matrix(0,2,n.step)
xy[,1] <- c(1/2,1/2)
```
```{r}
for(i in 2:n.step){
dx <- p*xy[1,i-1] + q*xy[1,i-1]*xy[2,i-1]
dy <- r*xy[1,i-1]*xy[2,i-1] + s*xy[2,i-1]
xy[,i] <- xy[,i-1] + c(dx,dy)*delta.t
}
plot(xy[1,],xy[2,],type="l",asp=TRUE)
```
平衡点では。
$$
\frac{\Delta x}{\Delta t} = p x + q xy = 0\\
\frac{\Delta y}{\Delta t} = r xy + s y = 0
$$
$$
x = -s/r
y = -p/q
$$
```{r}
xy <- matrix(0,2,n.step)
xy[,1] <- c(-s/r,-p/q)
for(i in 2:n.step){
dx <- p*xy[1,i-1] + q*xy[1,i-1]*xy[2,i-1]
dy <- r*xy[1,i-1]*xy[2,i-1] + s*xy[2,i-1]
xy[,i] <- xy[,i-1] + c(dx,dy)*delta.t
}
plot(xy[1,],xy[2,],asp=TRUE)
matplot(t(xy),type="l")
```
平衡点付近。
```{r}
xy <- matrix(0,2,n.step)
xy[,1] <- c(-s/r,-p/q)+ rnorm(2,0,0.01)
for(i in 2:n.step){
dx <- p*xy[1,i-1] + q*xy[1,i-1]*xy[2,i-1]
dy <- r*xy[1,i-1]*xy[2,i-1] + s*xy[2,i-1]
xy[,i] <- xy[,i-1] + c(dx,dy)*delta.t
}
plot(xy[1,],xy[2,],asp=TRUE)
matplot(t(xy),type="l")
```
まん丸に近い。
データから「まん丸の中心を計算してみる」。
```{r}
mean(xy[1,])
mean(xy[2,])
```
平衡点に近い。
```{r}
-s/r
-p/q
```
計算した「中心」を信じて微分(差分)方程式を、行列で表せるようにする。
中心を原点にする。
```{r}
new.x <- xy[1,] - mean(xy[1,])
new.y <- xy[2,] - mean(xy[2,])
new.xy <- rbind(new.x,new.y)
```
このnew.xyなら、その変化を表す行列が求められる。
```{r}
delta.xy12 <- (new.xy[,2:3] - new.xy[,1:2])
xy12 <- new.xy[,1:2]
xy12.inv <- solve(xy12)
(delta.xy12 %*% xy12.inv)/delta.t
```
```{r}
A <- 3:4
B <- 2:3
delta.xy12 <- (new.xy[,A] - new.xy[,B])
xy12 <- new.xy[,B]
xy12.inv <- solve(xy12)
(delta.xy12 %*% xy12.inv)/delta.t
```
```{r}
A <- 5:6
B <- A-1
delta.xy12 <- (new.xy[,A] - new.xy[,B])
xy12 <- new.xy[,B]
xy12.inv <- solve(xy12)
(delta.xy12 %*% xy12.inv)/delta.t
```
