---
title: CC2 modélisation
author: Oscar Plaisant
documentclass: scrartcl
header-includes: |
    \usepackage[top=2cm, bottom=2.5cm, left=2cm, right=2cm]{geometry}
    \usepackage{amsmath, amssymb, amsfonts, mathrsfs}
---

# Partie 1

## 1.

$z_1'$ est de classe $C^{1}$, donc on peut appliquer le théorème de Cauchy-Lipschitz.
On remarque que $z_1 = z_2 = 0 \mapsto 0$ est une solution du système.
Donc, par unicité de la solution, et par continuité de $z_1$, on sait que $z_1$ ne peut pas changer de signe, car sinon elle croiserait une autre solution.
On peut faire le même raisonnement pour $z_2$.
Comme on sait que $z_1^{0} > 0$ et $z_2^{0} > 0$, on sait que les deux fonctions seront toujours positives.


## 2.

$$
\begin{align}
(z_1 + z_2)' &= z_1' + z_2'\\
&= z_1(a_{12}z_2 - (a_{12}+a_{21})z_1z_2) + z_2(a_{21}z_2 - (a_{12} +a_{21})z_1z_2) \\
&= z_1z_2a_{12}-z_1z_2a_{12}-z_1z_2a_{21} +(z_2)^{2}a_{21} -z_1(z_2)^{2}a_{12}-z_1(z_2)^{2}a_{21} \\
&= -z_1z_2a_{21}+(z_2)^{2}a_{21}-(z_2)^{2}(a_{21}-(a_{12}+a_{21})z_1)
\end{align}
$$

## 3.

Comme on sait que $z_1+z_2 = 1$, et que les deux fonctions sont positives, on sait qu'elles sont minorées par 0, et si l'une devient plus grande que $1$ pour un $t$ donné, on sait que l'autre devrait être négative, ce qui est impossible, donc elles sont majorées par $1$.

On sait que les solutions existent tout le temps, car elles sont continues, et car elles sont bornées, donc leurs limites en $+\infty$ et $-\infty$ sont réelles.

## 4.

$z_1+z_2 = 1 \iff z_2 = 1-z_1$, donc on peut réécrire le système comme :  $z_1' = z_1(a_{12}(1-z_1)-(a_{12}+a_{21})z_1(1-z_1))$, et on a $z_2' = -z_1'$

Voici la version python de cette fonction :

```python
def z_prime(z: float, a_12: float, a_21: float) -> float:
    return z * (a_12 * (1-z) - (a_12 + a_21) * z * (1-z))
```

## 5.

Les points d'équilibre sont les points pour lesquels $z_1'(t) = 0$. $z_1' = z_1(a_{12}(1-z_1)-(a_{12}+a_{21})z_1(1-z_1))$, donc soit $z_1 = 0$, soit $a_{12}(1-z_1)-(a_{12}+a_{21})z_1(1-z_1) = 0$
$$
\begin{align}
\begin{cases}
z_1 = 0\\
a_{12}(1-z_1)-(a_{12}+a_{21})z_1(1-z_1) = 0 \\
\end{cases}
&\iff
\begin{cases}
z_1=0\\
z_1 = 1\\
a_{12}-z_1(a_{12}+a_{21}) = 0\\
\end{cases}\\
&\iff \begin{cases}
z_1 = 0 \\
z_1 = 1 \\
z_1 = \dfrac{a_{12}}{a_{12}+a_{21}}
\end{cases}
\end{align}
$$

Les points d'équilibre positifs ou nuls sont $t = 0$, $t = 1$ et $t = \dfrac{a_{12}}{a_{12} +a_{21}}$ si on suppose que cette valeur est dans $[0, 1]$.

## 6.

Pour la méthode d'Euler implicite, on doit résoudre cette équation :
$$
\begin{align}
z_1^{t+h} = z_1^{t}+h\cdot z_1'(z_1^{t+h}) &\iff z_1^{t+h} = z_1^{t}+h\cdot [ z_1^{t+h} (a_{12}(1-z_1^{t+h}) - (a_{12}+a_{21})z_1^{t+h}(1-z_1^{t+h})) ] \\
&\iff z_1^{t} + hz_1^{t+h}(1-z_1^{t+h}) [ a_{12} - z_1^{t+h}(a_{12}+a_{21}) ] - z_1^{t+h} = 0 \\
&\iff z_1^{t} + ha_{12}z_1^{t+h}(1-z_1^{t+h}) -h(a_{12}+a_{21})(z_1^{t+h})^{2}(1-z_1^{t+h}) - z_1^{t+h} = 0 \\
&\iff z_1^{t} + ha_{12}z_1^{t+h} -ha_{12}(z_1^{t+h})^{2} -h(a_{12}+a_{21})(z_1^{t+h})^{2} +h(a_{12}+a_{21})(z_1^{t+h})^{3} -z_1^{t+h} = 0 \\
&\iff h(a_{12}+a_{21})(z_1^{t+h})^{3} - h(2a_{12}+a_{21})(z_1^{t +h})^{2} + ha_{12}z_1^{t+h} + z_1^{t} = 0
\end{align}
$$

Comme on sait que $z_1^{t+h} \in [0, 1]$, on sait que l'on peut utiliser un développement limité comme approximation.
On cherche alors à résoudre :

$-h(2a_{12}+a_{21})(z_1^{t+h})^{2} +ha_{12}z_1^{t+h}+z_1^{t} = 0$
$\Delta = h^{2}a_{12}^{2} + 4z_1^{t}h(2a_{12}+a_{21})$

Les solutions de cette équation sont donc : 


$\displaystyle z_1^{t+h} = \frac{ha_{12} - \sqrt{ h^{2}a_{12}^{2}+4hz_1^{t}(2a_{12}+a_{21}) }}{2h(2a_{12}+a_{21})}$ et $\displaystyle z_1^{t+h} = \frac{ha_{12} + \sqrt{ h^{2}a_{12}^{2}+4hz_1^{t}(2a_{12}+a_{21}) }}{2h(2a_{12}+a_{21})}$

$\Delta > (ha_{12})^{2}$ donc $\sqrt{ \Delta } > ha_{12}$ par croissance de la racine carrée. On sait donc que $\displaystyle \frac{ha_{12}-\sqrt{ \Delta }}{2h(2a_{12}+a_{21})} < 0$ (on suppose $a_{12}$ et $a_{21}$ positifs).

Comme on cherche une solution positive, on sait que la solution à utiliser est : $\displaystyle z_1^{t+h} = \frac{ha_{12} + \sqrt{ h^{2}a_{12}^{2}+4hz_1^{t}(2a_{12}+a_{21}) }}{2h(2a_{12}+a_{21})}$

On peut donc coder une fonction qui calcule $z_1^{t+h}$ en fonction de $z_1^{t}$, $h$, $a_{12}$ et $a_{21}$ :

```python
def next_z_implicit(z: float, a_12: float, a_21: float, h: float) -> float:
    return (h*a_12 + np.sqrt(h**2*a_12**2 + 4*h*z*(W := 2*a_12 + a_21))) / (2*h*W)
```


Le code complet pour la simulation est donc :

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from typing import Callable

gradient = mpl.colormaps['plasma']


def z_prime(z: float, a_12: float, a_21: float) -> float:
    return z * (a_12 * (1-z) - (a_12 + a_21) * z * (1-z))

def next_z(z: float, a_12: float, a_21: float, h: float) -> float:
    return z + h * z_prime(z, a_12, a_21)

def next_z_implicit(z: float, a_12: float, a_21: float, h: float) -> float:
    return (h*a_12 + sqrt(h**2*a_12**2 + 4*h*z*(W := 2*a_12 + a_21))) / (2*h*W)

def euler_method(next_step: Callable, z_0, a_12, a_21, nb_steps: int, h=0.01):
    """The color changes depending on a_12 and a_21, but not z_0.
    This is so you can easily see what happens for any initial condition with
    the same parameters."""
    X = np.arange(nb_steps) * h
    values = []
    z = z_0
    for _ in X:
        values.append(z := next_z(z, h, a_12, a_21))

    plt.plot(values, linewidth=.6, color=gradient((a_12 + a_21) / 2))
    plt.xticks(range(nb_steps))


next_step_function = next_z
# next_step_function = next_z_implicit

# Plot for various initial conditions and values of a_12 and a_21
Z_0 = np.linspace(0, 1, 40)
for a_12 in np.linspace(.1, .9, 2):
    for a_21 in np.linspace(.1, .9, 5):
        euler_method(next_step_function,
                     z_0=Z_0, a_12=a_12, a_21=a_21,
                     nb_steps=4000, h=H)
plt.show()
```


## 7.

On observe que, sauf si les conditions initiales sont aux extrêmes ($z_1^{0} = 1$ ou $z_1^{0}=0$), les solutions convergent toutes vers un point d'équilibre. On peut donc conjecturer que $0$ et $1$ sont des points d'équilibre répulsifs, et que $\dfrac{a_{12}}{a_{12}+a_{21}}$ est le point d'équilibre attractif.

# Partie 2

## 8.

On obtient le système suivant :

$$\begin{cases}
z_1' = z_1(a_{12}z_2 + a_{13}(1-z_1-z_2) - q(z_1, z_2, 1-z_1-z_2)) \\
z_2' = z_2(a_{21}z_1 + a_{23}(1-z_1-z_2) - q(z_1, z_2, 1-z_1-z_2))
\end{cases}$$
## 9.

On commence par coder les fonctions correspondant à $z_1'$, $z_2'$ et $z_3'$ :

```python
def q(Z: list[float], A: mat_3x3) -> float:
    return (A[0][1]+A[1][0])*Z[0]*Z[1] + (A[0][2]+A[2][0])*Z[0]*Z[2] + (A[1][2]+A[2][1])*Z[1]*Z[2]

def z_1_prime(Z: list[float], A: mat_3x3) -> float:
    return Z[0] * ( A[0][1]*Z[1] + A[0][2]*Z[2] + q(Z, A) )

def z_2_prime(Z: list[float], A: mat_3x3) -> float:
    return Z[1] * ( A[1][0]*Z[0] + A[1][2]*Z[2] + q(Z, A) )

def z_3_prime(Z: list[float], A: mat_3x3) -> float:
    return Z[2] * ( A[2][0]*Z[0] + A[2][1]*Z[1] + q(Z, A) )
```

On code ensuite la fonction d'approximation par méthode d'Euler :

```python
def next_z(Z: list[float], A: mat_3x3, h: float) -> list[float]:
    return [Z[0] + h * z_1_prime(Z, A),
            Z[1] + h * z_2_prime(Z, A),
            Z[2] + h * z_3_prime(Z, A)]

def Replicator3G(A: mat_3x3, Z_0: list[float], t_max: float, h: float) -> list[float]:
    values = [Z_0]
    Z = Z_0
    for _ in np.arange(0, t_max, step=h):
        values.append(Z := next_z(Z, A, h))
    return values
```

Finalement, on code la fonction d'affichage de l'approximation :

```python
def plot_3G_replicator(A, Z_0, t_max: float, h: float) -> list[float]:
    values = Replicator3G(A, Z_0, t_max, h)
    values = np.transpose(np.array(values))
    col = gradient(np.mean(Z_0))  # just a number that changes depending on A
    plt.plot(values[0], linewidth=.6, color=col)
    plt.plot(values[1], linewidth=.6, color=col)
    plt.plot(values[2], linewidth=.6, color=col)
```

On obtient donc le code suivant :

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from typing import Callable, NewType

mat_3x3 = NewType('3x3 matrix', list[list[float]])

gradient = mpl.colormaps['plasma']

T = 3
h = 0.001
Z_1 = np.linspace(0, 1, 5)
Z_2 = np.linspace(0, 1, 5)
Z_3 = 1 - Z_1 - Z_2


def q(Z: list[float], A: mat_3x3) -> float:
    return (A[0][1]+A[1][0])*Z[0]*Z[1] + (A[0][2]+A[2][0])*Z[0]*Z[2] + (A[1][2]+A[2][1])*Z[1]*Z[2]

def z_1_prime(Z: list[float], A: mat_3x3) -> float:
    return Z[0] * ( A[0][1]*Z[1] + A[0][2]*Z[2] - q(Z, A) )

def z_2_prime(Z: list[float], A: mat_3x3) -> float:
    return Z[1] * ( A[1][0]*Z[0] + A[1][2]*Z[2] - q(Z, A) )

def z_3_prime(Z: list[float], A: mat_3x3) -> float:
    return Z[2] * ( A[2][0]*Z[0] + A[2][1]*Z[1] - q(Z, A) )

def next_z(Z: list[float], A: mat_3x3, h: float) -> list[float]:
    return [Z[0] + h * z_1_prime(Z, A),
            Z[1] + h * z_2_prime(Z, A),
            Z[2] + h * z_3_prime(Z, A)]

def Replicator3G(A: mat_3x3, Z_0: list[float], t_max: float, h: float) -> list[float]:
    values = [Z_0]
    Z = Z_0
    for _ in np.arange(0, t_max, step=h):
        values.append(Z := next_z(Z, A, h))
    return values


def plot_3G_replicator(A, Z_0, t_max: float, h: float) -> list[float]:
    values = Replicator3G(A, Z_0, t_max, h)
    values = np.transpose(np.array(values))
    col = gradient(np.mean(Z_0))  # just a number that changes depending on A
    plt.plot(values[0], linewidth=.6, color=col)
    plt.plot(values[1], linewidth=.6, color=col)
    plt.plot(values[2], linewidth=.6, color=col)

A_1 = [[0, 1, 1],
       [1, 0, 1],
       [1, 1, 0]]

A_2 = [[0, 1, 1],
       [1, 0, 1],
       [1, 1, 0]]

A_3 = [[ 0,  1, -1],
       [ 1,  0, -1],
       [-1, -1,  0]]

A_4 = [[ 0,  1, 1],
       [ 2,  0, 2],
       [-1, -1, 0]]

A_5 = np.array(A_4) * T


for z_1 in Z_1:
    for z_2 in Z_2:
        for z_3 in Z_3:
            plot_3G_replicator(
                A_5,  # change the parameters' matrix
                (z_1, z_2, z_3), T, h)

plt.show()
```