07-ODE_multidimension.Rmd

# ODE in multi-dimensions

## Terminology

$$\dot{X} = f(X, t, \Theta)$$

where:

- `X` is a vector of state variables (*e.g.* population size of several species)
- `t` stands for time
- `Θ` includes all the parameters
- `f` is a multi-dimensional function



## Linear system

Function `f` and its parameters `\theta`.
<!-- degeneracy -->






### General case

$$
\left\{
\begin{array}{l}
  \dot{x}_1 = \sum_i a_{1,i}x_i\\
  ... \\
  \dot{x}_n = \sum_i a_{n,i}x_{i} \\
  X(0) = X_0
\end{array}
\right.
$$

Using matrices:

$$\dot{X} = AX$$

Linear system in $n$ dimensions.

### Trivial case

If $i \neq j \Rightarrow 0$ then A in diagonal:


$$
\left\{
\begin{array}{l}
  \dot{x}_1 = a_{1,1}x_1\\
  ... \\
  \dot{x}_n = a_{1,1}x_{n} \\
  X(0) = X_0
\end{array}
\right.
$$

Then, all state variables are independent and the answer is obtained by solving
$n$ times a one-dimensional problem (described in the previous chapter).



### Diagonalizations

Now let's go back to the general case and assume that A is a diagonalizable
matrix (see Diagonalization on [wikipedia](https://en.wikipedia.org/wiki/Diagonalizable_matrix#Diagonalization)),
*i.e.* that it exits: $P$ and $D$ such as:

$$A = PDP^{-1}$$

This operation is a [change of basis](https://en.wikipedia.org/wiki/Change_of_basis),
we actually assume that it exists a basis in which our problem is as simple as the trivial
case, $P$ is the matrix that does the conversion. If our


### Graphical examples

```{r loadlinear, echo = F}
<<linear>>
```

```{r linear1, fig.width = 9, fig.height = 8}
beta <- list(
  matrix(c(-1, 0, 0, -1), 2),
  matrix(c(1, 0, 0, -1), 2),
  matrix(c(-1, 0, 0, 1), 2),
  matrix(c(1, 0, 0, 1), 2)
)
par(mfrow = c(2,2))
for (i in 1:4) intLinear(beta = beta[[i]])
```

```{r linear2, fig.width = 9, fig.height = 8}
beta <- list(
  matrix(c(0, -1, -1, 1), 2),
  matrix(c(0, 1, -1, 0), 2),
  matrix(c(0, -1, 1, 0), 2),
  matrix(c(0, 1, 1, 0), 2)
)
par(mfrow = c(2,2))
for (i in 1:4) intLinear(beta = beta[[i]])
```

```{r linear3, fig.width = 9, fig.height = 8}
beta <- list(
  matrix(c(0, 1, -1, -1), 2),
  matrix(c(0, 1, -1, 1), 2),
  matrix(c(1, 1, -1, 1), 2),
  matrix(c(-1, 1, -1, 1), 2)
)
par(mfrow = c(2,2))
for (i in 1:4) intLinear(beta = beta[[i]])
```

```{r linear4, fig.width = 9, fig.height = 8}
beta <- list(
  matrix(c(0, 1, -1, -2), 2),
  matrix(c(0, 1, -1, -1), 2),
  matrix(c(0, 1, -1, -.5), 2),
  matrix(c(0, 1, -1, -.1), 2)
)
par(mfrow = c(2,2))
for (i in 1:4) intLinear(beta = beta[[i]], mxt =100)
```

```{r linear5, fig.width = 9, fig.height = 8}
beta <- list(
  matrix(c(0, 2, -1, -.1), 2),
  matrix(c(0, 1, -1, -.1), 2),
  matrix(c(0, .5, -1, -.1), 2),
  matrix(c(0, .1, -1, -.1), 2)
)
par(mfrow = c(2,2))
for (i in 1:4) intLinear(beta = beta[[i]], mxt =100)
```



## Non-Linear systems


### Finding the eigen values

$$AX = \lambda X$$

Solutions of:

$$\left(A - \lambda I_n\right) X  = 0_n$$




## Examples

> Symetry drives weird shit! (K. S. McCann)


### Lotka Volterra

#### Intro

Classical model see @volterra_fluctuations_1927.

Equations:

$$
\left\{
\begin{array}{l}
  \dot{R} = rR - aCR \\
  \dot{C} = eaCR - mC
\end{array}
\right.
$$


#### First find equilibria and isolclines

$$
\left\{
\begin{array}{l}
  \dot{R} = 0 \Leftrightarrow C = r/a; C*=r/a \\
  \dot{C} = 0 \Leftrightarrow R = m/(ea); R* = m/(ea)
\end{array}
\right.
$$


#### Jacobian matrix

\begin{pmatrix}
r-aC* & -aR* \\
eaC* & eaR*-m
\end{pmatrix}

at the non-trivial equilibrium:

\begin{pmatrix}
0 & -aR* \\
eaC* & 0
\end{pmatrix}


#### Diagonalisation


\begin{pmatrix}
-\lambda & -aR* \\
eaC* & -\lambda
\end{pmatrix}

To get the two eigen values we must solve:

$$\lambda^2 + mr = 0$$


$$eig = \pm i \sqrt{mr}$$

(pure imaginary, so => spins)


#### Lotka-Volterra with`odeintr`

[`odeint`](http://headmyshoulder.github.io/odeint-v2/) is a C++ library, Timothy H. Keitt (may ring the bell to some of you)
has created a R package that wrapps around this library: [odeintr](https://github.com/thk686/odeintr).

Example from the GitHub project page:


```{r}
library(odeintr)
dxdt = function(x, t) c(x[1] - x[1] * x[2], x[1] * x[2] - x[2])
obs = function(x, t) c(Prey = x[1], Predator = x[2], Ratio = x[1] / x[2])
system.time({x = integrate_sys(dxdt, rep(2, 2), 20, 0.01, observer = obs)})
plot(x[, c(2, 3)], type = "l", lwd = 2, col = "steelblue", main = "Lotka-Volterra Phase Plot")
```



### Rozenzweig-McArthur model

#### Into

Another clasical, see @rosenzweig_graphical_1963 and @rosenzweig_paradox_1971
and McCann's book!

$$
\left\{
\begin{array}{l}
  \dot{R} = rR(1-\frac{R}{K}) - \frac{aCR}{R+R0} \\
  \dot{C} = \frac{eaCR}{R+R0} - mC
\end{array}
\right.
$$



<!-- 2Bedone
transcritical happens when one equilibrimu always there
Stable => unstabe
recepies for transcritical!


TC => unstable stable
Hopf Bofirctution =>  imaginary / cycle for ever
hop to K


Do the bifurcation diagramm K against Cmin Cmax -->