- 以下のテキストファイルはRmd文書。その整形はこちらを参照
- こちらにhtmlファイル
---
title: "ILAS2017"
author: "Ryo Yamada"
date: "2017年3月1日"
output:
html_document:
toc: true
toc_depth: 6
number_section: true
---
```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```
* コンピュータ言語処理は「コンピュータ上」で行う
* オープンソースと利用者・提供者コミュニティ
* 言葉は「まね〜コピー＆ペースト〜」が基本
* Windowsの場合: https://oku.edu.mie-u.ac.jp/~okumura/stat/R-win.html
* Macの場合: http://aoki2.si.gunma-u.ac.jp/R/begin.html
- アイコンをクリックして立ち上げる
- プロンプトに " 3 + 5 " と書き込み、"リターン・エンター"する
```{r}
3+5
```
- "8"が返ってくることを確認する
- 終了してみる。"q()" と入れて"リターン・エンター"する。qはquitの頭文字。()はおまじないのようなもの
```{r eval=FALSE}
q()
```
- 終了の仕方に２通りある
- 実行結果を保存する場合
- 実行結果を保存しない場合(このILASセミナーでは、基本的にはこちらで通します)
- 保存して終了してみる。
```{r}
x <- 3
y <- 5
z <- x + y
print(z)
```
```{r,eval=FALSE}
q()
```
- 再び起動
```{r}
print(x)
print(y)
print(z)
```
- 終了したときの状態から再開できる
- コマンドはすべてテキストファイルとして保存する
- テキストエディタ：フリーウェアでもよいのでひとつ選ぼう
- https://ferret-plus.com/810
```{}
x <- 3
y <- 5
z <- x + y
```
- 保存して、再利用するときは、コピーペースト(ファイルをRに読み込む方法もある)
- コマンドらしいところをコピーペーストしてみよう: http://d.hatena.ne.jp/ryamada/20170226
- コンピュータもコンピュータ言語も『道具』なので、自分がやりたいことができさえすればそれでよい。今、必要でないことはあえて覚えず、やりたいようにやりましょう。
- 計算
```{r}
3 + 5
3 * 5
3 - 5
3 / 5
3 ^ 5
```
- 変数に値を入れる
```{r}
x1 <- 3
x2 <- 5
```
- 入った値を確かめてみる
```{r}
x1
x2
```
- 変数を使って計算してみる
```{r}
x1 + x2
x1 - x2
x1 * x2
x1 / x2
x1 ^ x2
```
$$
y = a\sin{x} + b
$$
- 等間隔の数値の列を作る
- プロットする
```{r}
a <- 2
b <- 40
x <- seq(from=0,to=10,by=0.1)
y <- a * sin(x) + b
plot(x,y)
```
```{r}
plot(x,y,type="l")
```
- $y = x^2 - 3 x + 2 $ を$ -3 \le x \le 3$の範囲で描け。
- $y = (x - 2) \times (x - 1)$と計算しても同じである。このグラフを描け。
- 円 $(x-1)^2 + (y-2)^2 = 1 $ を描け。変数$t$を0から$2\pi$まで細かくとり、
$x-1 = \cos{t}$, $ y-2 = \sin{t}$ とするとよい。
```{r}
print(pi)
```
- 等間隔点を組み合わせることで格子状の点座標が作れる。
- 二次元格子を作る
```{r}
x <- seq(from = -3, to= 3 , by = 0.1)
y <- seq(from = -3, to = 3, by = 0.1)
# xとyとを組み合わせることで格子点のセットを作る
xy <- expand.grid(x,y)
head(xy)
X <- xy[,1]
Y <- xy[,2]
plot(X,Y,pch=20,cex=0.1)
```
- xy平面上の点に値を与え、その曲面を描く。
$$
z = x^2+y^2
$$
- 3次元プロットをするには、Rの機能を追加する。
- Rの追加機能はパッケージで行う
```{r,eval=FALSE}
install.packages("rgl")
library(rgl)
```
```{r, rgl = TRUE}
Z <- X^2+Y^2
plot3d(X,Y,Z)
```
2-1
$$
z = e^{-(x^2+y^2)}sin(x)*cos(y)
$$
なる曲面を描け。
```{r echo=FALSE}
x <- seq(from = -3, to= 3 , by = 0.1)
y <- seq(from = -3, to = 3, by = 0.1)
xy <- expand.grid(x,y)
head(xy)
X <- xy[,1]
Y <- xy[,2]
plot(X,Y)
Z <- exp(-(X^2+Y^2))*sin(X)*cos(Y)
plot3d(X,Y,Z)
```
$$
z = x^2+y^2
$$
- xとyの値を与えるとそれに応じたzという値を返す「関数」。
- コンピュータ言語でも、何かを与えて、それに応じて何かを返すものを関数と呼ぶ。
- 上式と同じ働きをする関数を作ってみる。
```{r,rgl=TRUE}
my.f1 <- function(x,y){
z <- x^2+y^2
return(z)
}
Z <- my.f1(X,Y)
plot3d(X,Y,Z)
```
- 与えるものを「引数」、返ってくるものを「返り値（かえりち）」と呼ぶ。
- 格子を作成するにあたり以下のようにした
```{r}
x <- seq(from = -3, to= 3 , by = 0.1)
y <- seq(from = -3, to = 3, by = 0.1)
# xとyとを組み合わせることで格子点のセットを作る
xy <- expand.grid(x,y)
```
- このように、from = a, to = b, by = kでxとyとをつくり、それを組み合わせた格子を返す関数は以下のように作れる。
```{r}
my.f2 <- function(a,b,k){
x <- seq(from = a, to = b, by =k)
y <- x
xy <- expand.grid(x,y)
return(xy)
}
A <- -4
B <- -2
K <- 0.5
gr <- my.f2(a=A,b=B,k=K)
plot(gr)
```
- ２次元平面の原点を中心とする半径rの円の円周上に等間隔に配意されたn個の点の座標を返す関数を作れ。ただし、等間隔点の発生には、等間隔角度を0から$2/pi$までつくり、それに対応する$\cos,\sin$を使って、$x,y$に対応する数値の列を作ればよい、数値列x,yを返すには、
```{r}
my.f3_1 <- function(r,n){
# x <- hoge
# y <- hoge
return(cbind(x,y)) #このようにすれば２列の行列で返る
# return(list(X=x,Y=y)) # このようにすれば、X,Yという２つのベクトルが返る
}
```
- 上で作成した関数( my.f3_1() )を呼び出して使う関数は次のようにして作れる。
次の関数my.f3_1_plot()は作成した座標をプロットする関数である。
```{r}
my.f3_1 <- function(r,n){
# hogehoge
}
my.f3_1_plot <- function(r,n){
tmp <- my.f3_1(r,n)
plot(tmp)
}
```
```{r echo = FALSE}
my.f3_1 <- function(r,n){
t <- seq(from=0,to=1,length=n) * 2*pi
x <- cos(t)
y <- sin(t)
return(cbind(x,y))
}
```
円周上に等間隔に配置した点(x,y)に$z = k \times sin(x\times y)$を与え、(x,y,z)をplot3d()関数を持ちいて３次元プロットする関数を作れ。
```{r echo=FALSE}
my.f3_2 <- function(k,r,n){
tmp <- my.f3_1(r,n)
z <- k * sin(tmp[,1]*tmp[,2])
plot3d(tmp[,1],tmp[,2],z)
}
my.f3_2(2,3,100)
```
$$
S_0 = 1\\
S_{k+1} = S_k + 1; k \ge 0
$$
$$
0,1,2,3,...
$$
- Rにこの数列を作らせる
```{r}
S0 <- 0
S1 <- S0 + 1
S2 <- S1 + 1
S0
S1
S2
```
- $S_0,...,S_n$の数列は$n+1$個の値の列
- Rでは値の列をベクトルという名前のオブジェクトで持つことができる
```{r}
S <- c(0,1,2)
print(S)
```
- 地道に計算して格納する
```{r}
S0 <- 0
S1 <- S0 + 1
S2 <- S1 + 1
S0
S1
S2
S <- c(S0,S1,S2)
S
```
- ベクトルの要素を確認する
- $S_2$は、ベクトルSの$2+1=3$番目の要素なので、次のように確認できる
```{r}
print(S[3])
```
- 数列を格納するベクトルを作る
```{r}
n <- 100
S <- rep(0,n+1)
print(S)
```
- 初期値設定
```{r}
S[1] <- 0
```
- 繰り返し・ループ
```{r}
for(i in 1:n){
S[i+1] <- S[i] + 1
}
print(S)
```
- 同じ仕事を何度もやりたいときには、『関数』を作る
```{r}
fsuuretsu1 <- function(n){
S <- rep(0,n+1)
S[1] <- 0
for(i in 1:n){
S[i+1] <- S[i] + 1
}
return(S)
}
fsuuretsu1
fsuuretsu1(5)
```
```{r}
fsuuretsu1(10)
```
- 数列 $S_{k+1} = S_{k} + 3; S_0 = 100$という数列の$S_0,S_1,S_2$を答えよ
- この数列の値$S_0,...,S_n$を長さn+1のベクトルで返す関数を作成せよ
$$
S_{k+2} = S_{k+1} + S_k\\
S_0 = 0\\
S_1 = 1
$$
- 地道な計算
```{r}
S0 <- 0
S1 <- 1
S2 <- S1 + S0
S3 <- S2 + S1
S4 <- S3 + S2
S5 <- S4 + S3
S <- c(S0,S1,S2,S3,S4,S5)
S
```
- ループを使う( for(i in 2:n) に注意)
```{r}
n <- 20
S <- rep(0,n+1)
S[1] <- 0
S[2] <- 1
for(i in 2:n){
S[i+1] <- S[i] + S[i-1]
}
print(S)
```
- フィボナッチ数$S_0,...,S_n$を長さn+1のベクトルで返す関数 ffibo() を作成せよ
```{r,echo=FALSE}
ffibo <- function(n){
S <- rep(0,n+1)
S[1] <- 0
S[2] <- 1
for(i in 2:n){
S[i+1] <- S[i] + S[i-1]
}
return(S)
}
```
$$
S_n = \frac{1}{\sqrt{5}}(\phi^n - (1-\phi)^n)\\
\phi = \frac{1+\sqrt{5}}{2}
$$
- 平方根
```{r}
phi <- (1+sqrt(5))/2
print(phi)
```
```{r}
n <- 4
1/sqrt(5)*(phi^n - (1-phi)^n)
```
- ベクトルを渡して、一気に計算する
```{r}
ns <- c(0,1,2,3,4,5,6,7,8,9,10)
1/sqrt(5)*(phi^ns - (1-phi)^ns)
```
- Rで整数列を作る
```{r}
n <- 10
ns <- 0:n
print(ns)
```
- 漸化式関数の結果と$\phi$を使った結果とを比較する
```{r}
n <- 20
S.zenka <- ffibo(n)
ns <- 0:n
S.phi <- 1/sqrt(5)*(phi^ns - (1-phi)^ns)
plot(S.zenka,S.phi)
abline(0,1,col=2)
```
$$
1 : \frac{1+\sqrt{5}}{2}
$$
- 近似解をRで探す
```{r}
x <- seq(from = -2, to = 2, length=10000)
y <- x^2 - x -1
plot(x,y,pch=20,cex=0.01)
abline(h=0,col=2)
```
- 近似解に対応するyの値は0に近いはず〜$|y-0|$がもっとも小さいxが近似解
- 絶対値、最小値、値を探す
```{r}
abs.y.0 <- abs(y-0)
plot(x,abs.y.0,pch=20,cex=0.1)
min(abs.y.0)
```
```{r}
dore <- which(abs.y.0==min(abs.y.0))
print(dore)
print(x[dore])
print((1+sqrt(5))/2)
plot(x,abs.y.0,pch=20,cex=0.1)
abline(v=x[dore],col=2)
```
$$
\phi = 1 + \cfrac{1}{1 + \cfrac{1}{1 + \cfrac{1}{1 + \cfrac{1}{\ddots}}}} = [1;1,1,1,1,\ldots]
$$
$$
A_1 = 1\\
A_2 = 1+\frac{1}{A_1}\\
A_3 = 1+\frac{1}{A_2}
$$
```{r}
n.iter <- 20
A <- rep(0,n.iter)
A[1] <- 1
for(i in 2:n.iter){
A[i] <- 1+1/A[i-1]
}
plot(A)
phi <- (1+sqrt(5))/2
abline(h=phi,col=2)
```
$$
\phi = \sqrt{1 + \sqrt{1 + \sqrt{1 + \sqrt{1 + \sqrt{\ldots}}}}}
$$
```{r}
n.iter <- 20
B <- rep(0,n.iter)
B[1] <- 1
for(i in 2:n.iter){
B[i] <- sqrt(1+B[i-1])
}
plot(B)
phi <- (1+sqrt(5))/2
abline(h=phi,col=2)
```
$$
\phi = 2 \cos{\frac{\pi}{5}}
$$
```{r}
2 * cos(pi/5)
phi
```
$$
e^{i\theta} + e^{-i\theta} = (\cos{\theta} + i \sin{\theta}) + (\cos{-\theta} + i \sin{-\theta})\\
= 2 \cos{\theta}
$$
であるから
$$
\phi = e^{i\frac{\pi}{5}} + e^{-i\frac{\pi}{5}} = 2 \cos{\frac{\pi}{5}}
$$
```{r}
print(1i)
print(1i * 1i)
exp(1i)
exp(1i * pi/5) + exp(-1i * pi/5)
```
- 複素表示を使ったフィボナッチ数
$$
\phi = e^{i\frac{\pi}{5}} + e^{-i\frac{\pi}{5}}\\
F(x) = \frac{1}{\sqrt{5}} (\phi^x - (1-\phi)^x)
$$
```{r, rgl = TRUE}
phi = exp(1i * pi/5) + exp(-1i * pi/5)
print(phi)
x <- seq(from=0,to=10,length=1000)
fx <- 1/sqrt(5) * (phi^x-(1-phi)^x)
plot(fx,type="l")
library(rgl)
plot3d(x,Re(fx),Im(fx),type="l")
fib <- ffibo(10)
spheres3d(0:10,fib,rep(0,length(fib)),radius=0.3,col=4)
```
- いわゆる、らせん
- 等角度の変化。
- 半径が一定速度で増加。
$$
x = r(t) \cos{\theta(t)}\\
y = r(t) \sin{\theta(t)}\\
t = 0,1,2,...\\
\theta(t) = a t\\
r(t) = b t
$$
```{r}
a <- pi/12
b <- 0.1
n <- 200
t <- 0:n
theta <- a * t
r <- b * t
x <- r*cos(theta)
y <- r*sin(theta)
plot(x,y,asp=TRUE,pch=20,type="b")
```
- 半径の伸びを緩やかにした、らせん
$$
r(t) = b \sqrt{t}
$$
```{r}
r <- b * sqrt(t)
x <- r*cos(theta)
y <- r*sin(theta)
plot(x,y,asp=TRUE,pch=20,type="b")
```
- 角度変化を黄金角にしてみる
$$
a = 2 \pi \times \frac{1+\sqrt{5}}{2}
$$
```{r}
a <- 2 * pi * (1 + sqrt(5))/2
theta <- a * t
r <- b * sqrt(t)
x <- r*cos(theta)
y <- r*sin(theta)
plot(x,y,asp=TRUE,pch=20)
```
- 領域に分けてみる(ボロノイ分割)
```{r}
# install.packages(c("deldir","sp")) # 二つのパッケージをインストールする
# http://stackoverflow.com/questions/9403660/how-to-create-thiessen-polygons-from-points-using-r-packages
# ボロノイ図を描く関数をウェブから借りてくる
voronoipolygons <- function(x) {
require(deldir)
require(sp)
if (.hasSlot(x, 'coords')) {
crds <- x@coords
} else crds <- x
z <- deldir(crds[,1], crds[,2])
w <- tile.list(z)
polys <- vector(mode='list', length=length(w))
for (i in seq(along=polys)) {
pcrds <- cbind(w[[i]]$x, w[[i]]$y)
pcrds <- rbind(pcrds, pcrds[1,])
polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
SP <- SpatialPolygons(polys)
voronoi <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
y=crds[,2], row.names=sapply(slot(SP, 'polygons'),
function(x) slot(x, 'ID'))))
}
X <- cbind(x,y) # ２次元座標
plot(X,pch=20,cex=0.1)
# ボロノイ多角形情報を取る
v <- voronoipolygons(X)
# ボロノイ多角形を描く
plot(v)
# 点を打つ
points(X,pch=20,cex=0.5,col="red") # 点の色はred
```
- 黄金角からずらしてみる
$$
a = 2 \pi \times (\frac{1+\sqrt{5}}{2} + \delta)
$$
```{r}
delta <- 0.5
a <- 2 * pi * ((1 + sqrt(5))/2 + delta)
theta <- a * t
r <- b * sqrt(t)
x <- r*cos(theta)
y <- r*sin(theta)
X <- cbind(x,y) # ２次元座標
plot(X,pch=20,cex=0.1)
# ボロノイ多角形情報を取る
v <- voronoipolygons(X)
# ボロノイ多角形を描く
plot(v)
# 点を打つ
points(X,pch=20,cex=0.5,col="red") # 点の色はred
```
$$
\frac{d^2 x(t)}{dt^2} = mg\\
x(0) = 0\\
v(0)=\frac{dx(0)}{dt} = 0
$$
- 微小時間後
$$
v(t+\Delta t) = v(t) + \frac{f(t)}{m} \times \Delta t = v(t) + g \Delta t\\
x(t+\Delta t) = x(t) + v(t) \times \Delta t
$$
$$
x(t) = \frac{1}{2}gt^2\\
v(t) = gt\\
f(t) = mg
$$
```{r}
maxt <- 5
deltat <- 0.5
t <- seq(from=0,to = maxt, by=deltat)
g <- 9.8
m <- 1
xt <- 1/2 * g * t^2
vt <- g * t
ft <- rep(m*g,length(t))
matplot(t,cbind(xt,vt,ft),type="l")
```
- 漸化式
$$
X_k = X_{k-1} + \Delta t \times V_k\\
V_k = V_{k-1} + \Delta t \times \frac{F_k}{m}\\
F_k = F_{k-1}\\
X_0 = 0\\
V_0 = 0\\
F_0 = mg
$$
```{r}
maxt <- 5
deltat <- 0.5
t <- seq(from=0,to = maxt, by=deltat)
X <- V <- F <- rep(0,length(t))
g <- 9.8
m <- 1
X[1] <- 0
V[1] <- 0
F[1] <- m*g
for(i in 2:length(t)){
X[i] <- X[i-1] + V[i-1] * deltat
V[i] <- V[i-1] + F[i-1]/m * deltat
F[i] <- F[i-1]
}
matplot(t,cbind(X,V,F),type="b",pch=20)
```
```{r}
matplot(cbind(xt,vt,ft,X,V,F),type="l")
```
$\Delta t$の値を変え、理論値と、漸化式の計算結果との違いを比較せよ。
$$
\frac{dx}{dt} = a x\\
x(0) = b
$$
```{r}
maxt <- 5
deltat <- 0.5
t <- seq(from=0,to = maxt, by=deltat)
x <- rep(0,length(t))
a <- 0.5
b <- 10
x[1] <- b
for(i in 2:length(t)){
x[i] <- x[i-1] + a * x[i-1] * deltat
}
plot(t,x,type="b",pch=20)
```
```{r}
maxt <- 5
deltat <- 0.001
t <- seq(from=0,to = maxt, by=deltat)
x <- rep(0,length(t))
a <- 0.5
b <- 10
x[1] <- b
for(i in 2:length(t)){
x[i] <- x[i-1] + a * x[i-1] * deltat
}
plot(t,x,type="b",pch=20)
```
$$
\frac{dx}{dt} = px + qy\\
\frac{dy}{dt} = rx + sy\\
x(0) = 1\\
y(0) = 1
$$
$$
\begin{pmatrix} \frac{dx}{dt}\\\frac{dy}{dt} \end{pmatrix} = \begin{pmatrix}p,q\\r,s \end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix}
$$
```{r}
maxt <- 5
deltat <- 0.00005
t <- seq(from=0,to = maxt, by=deltat)
x <- y <- rep(0,length(t))
x[1] <- 1
y[1] <- 1
p <- 1
q <- -2
r <- 3
s <- -1
for(i in 2:length(t)){
x[i] <- x[i-1] + (p * x[i-1] + q * y[i-1]) * deltat
y[i] <- y[i-1] + (r * x[i-1] + s * y[i-1]) * deltat
}
par(mfcol=c(1,2))
matplot(t,cbind(x,y),type="l")
plot(x,y,type="l")
eigen(matrix(c(p,q,r,s),2,2))
```
```{r}
maxt <- 5
deltat <- 0.00005
t <- seq(from=0,to = maxt, by=deltat)
x <- y <- rep(0,length(t))
x[1] <- 1
y[1] <- 1
p <- 1
q <- -2
r <- 3
s <- -1.5
for(i in 2:length(t)){
x[i] <- x[i-1] + (p * x[i-1] + q * y[i-1]) * deltat
y[i] <- y[i-1] + (r * x[i-1] + s * y[i-1]) * deltat
}
par(mfcol=c(1,2))
matplot(t,cbind(x,y),type="l")
plot(x,y,type="l")
eigen(matrix(c(p,q,r,s),2,2))
```
- p,q,r,sの値をいろいろ変えて、プロットしてみよ
- 変数の数を３つに増やして、挙動を調べてみよ
- 被捕食動物の増減は、自然増の項と被捕食による減少の項からなる
- 捕食動物の増減は、自然減の項と捕食による増加の項からなる
$$
\frac{dx}{dt} = p x - q xy\\
\frac{dy}{dt} = r xy - s y\\
p,q,r,s > 0
$$
```{r}
maxt <- 20
deltat <- 0.001
t <- seq(from=0,to = maxt, by=deltat)
x <- y <- rep(0,length(t))
x[1] <- 2
y[1] <- 1
p <- 1
q <- 1
r <- 1
s <- 1
for(i in 2:length(t)){
x[i] <- x[i-1] + (p * x[i-1] - q * x[i-1] * y[i-1]) * deltat
y[i] <- y[i-1] + (r * x[i-1] * y[i-1] - s * y[i-1]) * deltat
}
par(mfcol=c(1,2))
matplot(t,cbind(x,y),type="l")
plot(x,y,type="l")
```
- 増減の無い状態
$$
\frac{dx}{dt} = p x - q xy = 0\\
\frac{dy}{dt} = r xy - s y = 0
$$
- 上記を満足するx,yでは、xもyも増減しない。そのことを、そのような初期値から開始することで確かめよ
- 増減しない(x,y)の値から少しだけずれた状態からのシミュレーション結果と、大幅にずれた状態からのシミュレーション結果を比較せよ
- 1次元空間の拡散
- n 個の小部屋が並んでいる
- 時刻tのとき、i番目の部屋には$x_i(t)$の「量」がある
- 単位時間$\Delta t$後に、$x_i(t)$のうち、$f_L\times x_i(t)$が左の小部屋に、$f_R \times x_i(t) が右の小部屋に移動し、$(1-f_L-f_R)\times x_i(t)$は同じ小部屋にとどまる
$$
x_i(t+\Delta t) = x_i(t) - f_L x_i(t) - f_R x_i(t) + f_L x_{i+1}(t) + f_R x_{i-1} (t)
$$
```{r}
n.room <- 10
x <- 1:n.room
maxt <- 500
deltat <- 1
t <- seq(from=0,to = maxt, by=deltat)
fL <- fR <- 0.01
X <- matrix(0,length(t),n.room)
X[1,] <- rnorm(n.room)
for(i in 2:length(t)){
for(j in 1:n.room){
out <- X[i-1,j] * (fL+fR)
in1 <- in2 <- 0
if(j>1){
in1 <- X[i-1,j-1] * fR
}
if(j<n.room){
in2 <- X[i-1,j+1] * fL
}
in.out <- in1 + in2 - out
X[i,j] <- X[i-1,j] + in.out
}
}
matplot(X,type="l")
```
- より広い１次元空間、矩形初期値
- 以下のコードをコピーペーストで実行してみよ
```{r}
# 部屋の数
n <- 1000
# 追跡時刻数
n.t <- 1000
# 砂の量を記録する行列
X <- matrix(0,n.t,n)
# 初期時刻の砂の量(高さ)をところどころに与える
# m箇所に砂がある
m <- n/10
# 砂のある場所
X[1,10:50] <- 10
X[1,200:300] <- 100
X[1,500:800] <- 50
# 隣の部屋へは、delta.t * D の割合で移動する
delta.t <- 1
D <- 0.4 # D は0.5以下
for(i in 2:n.t){
# t = iの部屋の砂の量を仮に決める
tmp <- X[i-1,]
# 全部の部屋の砂の移動を、一部屋ずつ処理する
for(j in 1:n){
to.Left <- delta.t * D * X[i-1,j]
to.Right <- delta.t * D * X[i-1,j]
# 自分の部屋の砂は減る
tmp[j] <- tmp[j] - to.Left - to.Right
# 左隣の部屋は増える
if(j-1 >= 1){
tmp[j-1] <- tmp[j-1] + to.Left
}
# 右隣の部屋は増える
if(j+1 <= n){
tmp[j+1] <- tmp[j+1] + to.Right
}
}
X[i,] <- tmp
}
image(X)
```
```{r eval=FALSE}
#par(ask=TRUE)
for(i in 1:n.t){
plot(1:n,X[i,],ylim = c(0,max(X)),type="l")
}
```
- 初期状態として、サインカーブ状のものを与え、拡散の様子を確認せよ
- 右向きと左向きとの移動量の割合を変えて実行してみよ
- 上下左右の小部屋に移動する
$$
x_{i,j}(t+\Delta t) = x_{i,j}(t) - a + b\\
a = f_L x_i(t) f_R x_i(t) + f_U x_i(t) + f_L x_i(t)\\
b = f_L x_{i+1,j}(t) + f_R x_{i-1,j} (t) + f_U x_{i,j-1} + f_B x_{i,j+1}
$$
```{r}
n.room <- 5
xy <- as.matrix(expand.grid(1:n.room,1:n.room))
plot(xy)
```
```{r}
n.room <- 40
maxt <- 50
deltat <- 1
t <- seq(from=0,to = maxt, by=deltat)
fL <- fR <- 0.1
fU <- 0.2
fB <- 0
X <- array(0,c(length(t),n.room,n.room))
xy <- as.matrix(expand.grid(1:n.room,1:n.room))
X[1,18:20,10:15] <- 1
X[1,30:35,20:25] <- 2
for(i in 2:length(t)){
for(j in 1:n.room){
for(k in 1:n.room){
out <- X[i-1,j,k] * (fL+fR+fU+fL)
in1 <- in2 <- in3 <- in4 <- 0
if(j>1){
in1 <- X[i-1,j-1,k] * fR
}
if(j<n.room){
in2 <- X[i-1,j+1,k] * fL
}
if(k>1){
in3 <- X[i-1,j,k-1] * fU
}
if(k<n.room){
in4 <- X[i-1,j,k+1] * fB
}
in.out <- in1 + in2 +in3 + in4 - out
X[i,j,k] <- X[i-1,j,k] + in.out
}
}
}
```
```{r}
for(i in c(1,5,10,15,20,25,30,35)){
image(X[i,,])
}
```
```{r,eval=FALSE}
for(i in 1:30){
#for(i in 1:length(t)){
image(X[i,,])
}
```
```{r,eval=FALSE}
a <- 2.8 * 10^(-4)
b <- 5 * 10^(-3)
tau <- 0.1
k <- -0.005
size <- 100
dx <- 2/size
T <- 10.0
dt <- 0.9 * dx^2/2
n <- floor(T/dt)
U <- matrix(runif(size^2),size,size)
V <- matrix(runif(size^2),size,size)
ctr <- 2:(size-1)
plus <- 3:size
minus <- 1:(size-2)
U.hx <- list()
V.hx <- list()
U.hx[[1]] <- U
V.hx[[1]] <- V
for(i in 1:n){
deltaU <- (U[minus,ctr] + U[plus,ctr] + U[ctr,minus] + U[ctr,plus] - 4 * U[ctr,ctr])/(dx^2)
deltaV <- (V[minus,ctr] + V[plus,ctr] + V[ctr,minus] + V[ctr,plus] - 4 * V[ctr,ctr])/(dx^2)
Uctr <- U[ctr,ctr]
Vctr <- V[ctr,ctr]
new.U <- U
new.V <- V
new.U[ctr,ctr] <- Uctr + dt * (a*deltaU + Uctr - Uctr^3 -Vctr + k)
new.V[ctr,ctr] <- Vctr + dt * (b*deltaV + Uctr - Vctr)/tau
new.U[1,] <- new.U[2,]
new.U[size,] <- new.U[size-1,]
new.U[,1] <- new.U[,2]
new.U[,size] <- new.U[,size-1]
new.V[1,] <- new.V[2,]
new.V[size,] <- new.V[size-1,]
new.V[,1] <- new.V[,2]
new.V[,size] <- new.V[,size-1]
U <- new.U
V <- new.V
if(i%%10000 == 1){
image(U,col=topo.colors(100))
U.hx[[i+1]] <- U
V.hx[[i+1]] <- V
}
}
```
```{r,eval=FALSE}
par(ask=TRUE)
for(i in 1:length(U.hx)){
if(is.matrix(U.hx[[i]])){
image(U.hx[[i]],col=topo.colors(100))
}
}
for(i in 1:length(V.hx)){
if(is.matrix(V.hx[[i]])){
image(V.hx[[i]],col=topo.colors(100))
}
}
par(ask=FALSE)
```
