---
title: "Résolution d'EDO avec une librairie"
author: "Thomas Lepoutre"
date: "2025-09-06"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

La résolution d'équations différentielles ordinaires (par opposition à EDP, partielles) est un problème important du calcul scientifique. De façon logique, des outils ont été développés pour ne pas refaire les schémas lorsque ce n'est pas nécessaire. L'outil qui nous intéresse ici est donc la librairie desolve. Elle s'utilise pour résoudre les problèmes de Cauchy $$
\begin{cases}
X'(t)=F(t,X(t))\\
X(t_{0})=X_0
\end{cases}
$$

Dans cet état dR'esprit $X:[t_0,t_max]\mapsto \mathbb{R}^d$ et \$F:]t_0,t\_{max}\times \mathbb{R}^d^\mapsto\mathbb{R}d \$ On considérera toujours dans la programmation que la fonction dépend de $t$ même si ce n'est pas le cas. On fait cela car les routines ne sont pas faites pour distinguer les deux cas. On regarde un exemple simple $$
\begin{cases}
x'(t)=-x\\
x(0)=1
\end{cases}
$$ Dont vous connaissez la solution $x(t)=e^{-t}$. Il faut commencer par faire appel à la libairie

```{r}
library(deSolve)
Linearneg <- function(t, X, parameters) {
  {
    dX <-  -2*X
    list(c(dX))
  }
}
```

On crée une liste de temps auxquels on veut calculer la fonction (ELLE DOIT CONTENIR $t_0$)

```{r}
t0=0;
x0=1;
X=x0;
times=t0+0.01*(0:200)
out <- ode(y = X, times = times, func = Linearneg)
plot(out)
```

Dans le même esprit, essayer de programmer la résolution du problème suivant

$$
\begin{cases}
x'(t)=x*(1-x)\\
x(0)=1/2
\end{cases}
$$

```{r}
# A vous 

```

Un peu plus dur: la simulation de système

$$
\begin{cases}
x'(t)=2x-xy,\\
y'(t)=xy-y,\\
x(0)=y(0)=1
\end{cases}
$$

```{r}
Lotka <- function(t, state, parameters) {
  {
    x <-  state[1]
    y <- state[2]
    dxdt<-2*x-x*y
    dydt<-x*y-y
    dX<-c(dxdt,dydt)
    list(dX)
  }
}
```

```{r}
t0=0;
x0=c(1,1);
X=x0;
times=t0+0.01*(0:1000)
out <- ode(y = X, times = times, func = Lotka)
plot(out)
plot(out[,2],out[,3])
plot(out[,1],out[,2])
par(new=TRUE)
plot(out[,1],out[,3])
```

ATTENTION: out[,1]est le temps (les autres variables suivent).

Une remarque sur ce modèle, les solutions sont périodiques. Est ce que ça peut se voir sur le deuxième graphe?

Pour continuer, il faut simuler un système SIR

$$
\begin{cases}
S'(t)=-SI,\\
I'(t)=SI-\mu I -\beta I,\\
R'(t)=\beta I
S(0)=10, I(0)=0.1, R(0)=0,\\
\mu=0.1, \beta=1
\end{cases}
$$

```{r}
# A vous
```

Le schéma d'Euler

Si on avait pas une routine desolve, il faudrait faire un schéma numérique. Le principe du schéma d'Euler est le suivant. On APPROCHE la solution de

$$
X'(t)=F(t,X(t))
$$

Aux temps $t_0, t_0+h, \dots t_0+nh$

par la suite

$$
Y_0 =X(t_0), Y_{k+1}=Y_k+ F(t_k, Y_k)*h. 
$$

C'est basé sur l'idée que

$$
X'(t)\sim \frac{X(t+h)-X(t)}{h}
$$

Et donc qu'en un certain sens, l'équation peut se voir comme

$$
\frac{X(t+h)-X(t)}{h}\sim F(t,X(t)). 
$$

C'est dans le TP boucles.