---
title: "ローレンツ方程式"
author: "宮部賢志"
date: "2015年9月29日"
output:
  html_document:
    number_sections: yes
    self_contained: no
    toc: yes
---

<link rel="stylesheet" href="https://rawgithub.com/brentertz/scianimator/master/assets/css/scianimator.css"/>
<script src="https://ajax.googleapis.com/ajax/libs/jquery/1.7.1/jquery.min.js"></script>
<script src="https://rawgithub.com/brentertz/scianimator/master/assets/js/jquery.scianimator.min.js"></script>

```{r setup, include=FALSE}
# アニメーション作成用のフック関数の指定
library(knitr)
opts_knit$set(animation.fun = hook_scianimator)
```

#ローレンツ方程式
*ローレンツ方程式*とは，カオス的振る舞いを示す非線形方程式の1つである．
$$\begin{align*}
\dfrac{dx}{dt}=&-px+py\\
\dfrac{dy}{dt}=&-xz+rx-y\\
\dfrac{dz}{dt}=&xy-bz
\end{align*}$$
ローレンツ(Lorentz)の1963年の論文の中で，
$p=10, r=28, b=8/3$という設定が紹介されている．

#シミュレーション
以下ではこのローレンツ方程式のシミュレーションを行う．
単純のためライブラリを使うが，レポートではできる限りライブラリを使わず，アルゴリズムの勉強をしてほしい．

必要なライブラリを読み込む．

```{r}
library(deSolve)
library(rgl)
```

パラメーターを設定する．
```{r}
parameters <- c(s=10, b=8/3, r=28)
initial <- c(x=0, y=0.1, z=0)
times <- seq(0, 20, 0.1)
```

ローレンツ方程式を解いて，その結果をoutに出力する．
```{r}
Lorentz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx <- s * (y-x)
    dy <- r * x - y - x * z
    dz <- -b * z + x * y
    list(c(dx,dy,dz))
  })
}
out <- ode(y = initial, times = times, func = Lorentz, parms = parameters)
```

その結果をcx,cy,czにコピーしながら動画として再生する．
```{r, lorentz-fig, fig.show="animate", interval=0.1, fig.width=4, fig.height=4}
imax <- nrow(out)
cx <- rep(NA,imax)
cy <- rep(NA,imax)
cz <- rep(NA,imax)

library(scatterplot3d)
xlim <- c(min(out[,"x"]),max(out[,"x"]))
ylim <- c(min(out[,"y"]),max(out[,"y"]))
zlim <- c(min(out[,"z"]),max(out[,"z"]))
for(i in 1:imax){
    cx[i] <- out[i,"x"]
    cy[i] <- out[i,"y"]
    cz[i] <- out[i,"z"]
    scatterplot3d(cx,cy,cz,xlim=xlim,ylim=ylim,zlim=zlim,type="l")
}
```
