#### Hugo CRESPO, Mélanie FLESSEL, Xinyun LI
# TD Intégration géométrique : Méthode d'Euler

L’objectif de ce TD est de comparer les formes classiques et les formes symplectiques de la méthode la plus élémentaire d’intégration numérique des EDO, le schéma d’intégration d’Euler, par des simulations de différents exemples de systèmes hamiltoniens.

Nous considérerons l’espace d’état $\mathsf{X} = \mathbb{R} \ni x$, un intervalle de temps $[0;T]$ $\ni$ t avec $T > 0$ et l’EDO d’ordre 1 :

$$
\frac{dx}{dt} = f(x)
$$
où $f$ est une fonction $C^{\infty}$

Nous noterons la condition initiale $x(0) = x_0$ dans la suite.
Les méthodes d’Euler font partie des méthodes “à un pas” et nous noterons $h$ le pas d’intégration, fixe. Par conséquent nous cherchons une solution approchée $x_k$ de la solution $x(k,h)$ de l'EDO au temps $(k,h)$ et choisissons $T=Nh$ avec $N$ un nombre entier.

## Méthode d'Euler explicite
La méthode d’Euler explicite a pour équation de récursion
$$
x_{k+1} = x_k + hf(x_k)
$$

### Implémentation d'Euler Explicit

function E=eulerExplicit(f,a,b,y0,M)
h=(b-a)/M;    % pas de temps
Y=zeros(2,M+1);  % creation d'une matrice de dimensions 2 x (M+1)
%T=a:h:b;
Y(1,1)=y0(1); % initialisation de p
Y(2,1)=y0(2); % innitilisation de q
for j=1:M
    r = f(Y(1,j),Y(2,j));     % appel de la fonction
    Y(1,j+1)=Y(1,j) + h*r(1); % systeme hamiltonien
    Y(2,j+1)=Y(2,j) + h*r(2);
end
plot(Y(1,:),Y(2,:))  % plot q par rapport à p
E = Y';
end

## Méthode d'Euler implicite
La méthode d’Euler implicite a pour équation de récursion
$$
x_{k+1} = x_k + hf(x_{k+1})
$$
### Implémentation d'Euler Implicite

% Backward Euler
function E=eulerImplicit(f,a,b,y0,M)
t = 0:1000; %time vector
dt = t(2)-t(1); %time step
Y(1,1)=y0(1);  % initialisation de p
Y(2,1)=y0(2);  % innitilisation de q
options = optimset('TolX' , 1e-6);  % option d'arret : 1e-6
for i = 1:length(t)-1  % resolution du systeme f_n+1 = (p_n+1,q_n+1)
    Y(:, i+1) = fsolve(@(y) Y(:, i) + h*f(y(1), y(2)) - y, Y(:, i),options);
end 
plot(Y(1,:),Y(2,:))  % plot q par rapport à p
E = Y';
end

## Méthode d'Euler symplectique : implicite en q, explicite en p
Adaptation des deux méthodes précédentes, cette méthode d’intégration consiste, en l’application de la méthode d’Euler implicite pour le calcul du déplacement généralisé $q$ et de la méthode d’Euler explicite pour le calcul du moment généralisée $p$.

### Implémentation d'Euler Symplectique

function E = eulerSymplectique(f,a,b,y0,M)
h=(b-a)/M;
t = 0:1000; %time vector
dt = t(2)-t(1); %time step
Y(1,1)=y0(1);  % initialisation de p
Y(2,1)=y0(2);  % innitilisation de q
options = optimset('TolX' , 1e-6);  % option d'arret : 1e-6
for i = 1:length(t)-1  % resolution du systeme f_n+1
    r = f(Y(1,i),Y(2,i));     % appel de la fonction
    Y(1,i+1)=Y(1,i) + h*r(1); % p_n+1
    % calcul de q_n+1
    Y(:, i+1) = fsolve(@(y) Y(:, i) + h*f(Y(1,i), y(2)) - y, Y(:, i),options);
end 
plot(Y(1,:),Y(2,:))  % plot q par rapport à p
E = Y';
end

## Tests

### Oscillateur Harmonique
L'oscillateur harmonique est défini par un système hamiltonien standard engendré par la fonction hamiltonienne :
$$
H_0(q, p) = K(p) + U(q)
$$
où $q$ est le déplacement généralisé et $p$ est la quantité de mouvement généralisée.

Pour l’application numérique, nous considèrons un système masse-ressort où
- $q$ est l’élongation du ressort, initialisée à 1
- $p$ est la quantité de mouvement, intialisée à 0
- $U(q)=\frac{1}{2}kq^2$, est l’énergie potentielle élastique du ressort
- $K(p)=\frac{1}{2}\frac{p^2}{m}$, est l’énergie cinétique de la masse
- $k = 214$ N/m, est la raideur du ressort
- $m$ = 4.7$ kg, la masse considérée

#### Euler explicit

f=@(p,q) [2*k*q ; -2*p/m];
a=0;
b=10;
ya = [0;1] % initialisation p=0, q=1
M=1000;
YY=eulerExplicit(f,a,b,ya,M) %YY=eulerImplicit(f,a,b,ya,M) 


#### Euler implicit

Par la méthode d'Euler explicite, on obtient une oscillation dont l'amplitude augmente, au lieu du mouvement périodique de la solution exacte.
A l'inverse, la méthode d'Euler implicite nous donne une oscillation dont l'amplitude diminue.
Ces méthodes sont instables pour l'oscillateur harmonique.

#### Euler symplectique
- 1000 pas de temps, entre 0 et 20
- Avec initialisation $p = 0$, $q = 1$

Par la méthode d'Euler symplectique, on obtient cette fois une oscillation dont l'amplitude converge et diverge, autour de la solution
Augmenter le pas, de temps nous permet finalement d'approcher la solution cherchée.

### Pendule simple
Le pendule simple en mouvement dans un plan est modélisé par le système hamiltonien :

$$
H_0(q, p) = K(p) + U(q)
$$

où $q$ est l’angle de la barre par rapport à la position de repos (verticale) et $p$ est la quantité de mouvement de la barre en rotation.

De plus, on considère les paramètres suivants :
- $U(q)=-mgl_0\cos(q)$, l'énergie potentielle de gravité
- $K(p)=\frac{1}{2} \frac{p^2}{J}$, l’énergie cinétique en rotation de la barre
- $g=9.81$ $m.s^{−2}$, la constante de gravité
- $m = 0.11$ kg, la masse de la barre
- $l_0 = 0, 15$ m, la position de son centre de gravité
- $J = 0, 24$ $kg.m^2$, son inertie de rotation

#### Euler explicit
- 1000 pas de temps, entre 0 et 10
- Avec initialisation $p = 0$, $q = 1$

J = 0.24
m = 0.11
l0 = 0.15
g = 9.8
f = @(p,q) [0 1; -1 0]*[p/J; m*g*l0*sin(q)]
a=0;
b=10;
ya = [0;1]
YY=eulerExplicit(f,a,b,ya,M) %YY=eulerImplicit(f,a,b,ya,M) 


- 1000 pas de temps, entre 0 et 40 (pour avoir plus de période)

- D'après le graphe qu'on a obtenu, on voit bien que la méthode **Euler Explicit** diverge, **s'ésloigne** vers extérieur.

#### Euler implicit
On utilise les mêmes condtions initiales que pour la méthode explicite.

- 1000 pas de temps, entre 0 et 10
- Avec initialisation $p = 0$, $q = 1$

% Initilisation des variables
J = 0.24
m = 0.11
l0 = 0.15
g = 9.8
f = @(p,q) [0 1; -1 0]*[p/J; m*g*l0*sin(q)]
a=0;
b=10;
ya = [0;1]
YY=eulerImplicit(f,a,b,ya,M)


- 1000 pas de temps, entre 0 et 20

- 1000 pas de temps, entre 0 et 40(plus de périodes)

- D'après les graphes qu'on a obtenus, on voit bien que la méthode **Euler Implicite** converge vers intérieur.

#### Euler symplectique
- 1000 pas de temps, entre 0 et 10
- Avec initialisation $p = 0$, $q = 1$

% Initilisation des variables
J = 0.24
m = 0.11
l0 = 0.15
g = 9.8
f = @(p,q) [0 1; -1 0]*[p/J; m*g*l0*sin(q)]
a=0;
b=10;
ya = [0;1]
M = 1000
YY=eulerSymplectique(f,a,b,ya,M) % appèle la fonction eulerSympletique


- 1000 pas de temps, entre 0 et 20

Finalement, pour de petits pas de temps, cette méthode est **stable** pour l'oscillateur harmonique et le pendule simple

## Lotka-Volterra
C’est un système hamiltonien d’ordre 2 définis par rapport un crochet de Poisson **non canonique** : il est de rang 2 sauf sur une sous-variété de codimension 1.
Nous allons considérer ici un modèle appelé proie-prédateur, qui modélise des sys-
tèmes écologiques.

C’est le système hamiltonien défini sur un ensemble ouvert de $\mathbb{R^2}$, l’espace d’état $\mathsf{X} = \mathbb{R^*_+} \times \mathbb{R^*_+}$ défini par la matrice de Poisson

$$
J(x)=
\begin{pmatrix}
0 & x_1x_2 \\
-x_1x_2 & 0
\end{pmatrix}
$$

et généré par la fonction hamiltonienne
$$
H_0(x) = x_1 − \ln(x_1) + x_2 − 2\ln(x_2)
$$

Sur matlab, on utilise dans ce cas la fonction :

f = @(x1,x2) (x1*x2)*[0 1;-1 0]*[1 - 1/x1; 1 - 2/x2]

Toutefois, dans le domaine $\mathsf{X}$, on peut définir les variables canoniques à partir des coordonnées naturelles par le changement de variables d’état suivant :

$$
\begin{pmatrix}
q \\
p
\end{pmatrix} = \begin{pmatrix}
\ln(x_1)\\
\ln(x_2)
\end{pmatrix}
$$

### Implémentation
#### Euler Explicit
- 1000 pas de temps, entre 0 et 10
- Avec initialisation $x_1 = 2$, $x_2 = 1$

#### Euler Implicit
- 1000 pas de temps, entre 0 et 10
- Avec initialisation $x_1 = 2$, $x_2 = 1$

#### Euler Symplectique
- 1000 pas de temps, entre 0 et 10
- Avec initialisation $x_1 = 2$, $x_2 = 1$

D'après les cours, on sait que la méthode de Euler Symplectique ne peut pas marcher avec la modèle de Lotka-Volterra, mais dans notre cas, d'après le graphe obtenu, la méthode marche.
En effet, dans notre fonction eulerSympletique, on a utilisé directement $(x_1x_2)* J_s*H_0(q, p)$ lors de l'appel de la fonction, ce qui a corrigé la méthode originale.

#### Lotka-Volterra en coordonnées naturelles

function E = eulerSymplectique(f,a,b,y0,M)
h=(b-a)/M;
t = 0:1000; %time vector
dt = t(2)-t(1); %time step
tic  
Y(1,1)=y0(1);  % initialisation de p
Y(2,1)=y0(2);  % innitilisation de q
options = optimset('TolX' , 1e-6);  % option d'arret : 1e-6
for i = 1:length(t)-1  % resolution du systeme f_n+1
    r = f(Y(1,i),Y(2,i));     % appel de la fonction
    Y(1,i+1)=Y(1,i) + h*r(1); % p_n+1
    Y(:, i+1) = fsolve(@(y) Y(:, i) + h*f(Y(1,i), y(2)) - y, Y(:, i),options);
end 
plot(Y(1,:),Y(2,:))  % plot 
hold on
plot(Y(1,:))  % plot de x1
plot(Y(2,:))  % plot de x2
hold off
E = Y';
end

 %appele fonction Lotka-Volterra original avec point de initialisation (x1,x2) = (2,5)
f = @(x1,x2) (x1*x2)*[0 1;-1 0]*[1 - 1/x1; 1 - 2/x2]
 appele fonction Lotka-Volterra original
f = @(x1,x2) (x1*x2)*[0 1;-1 0]*[1 - 1/x1; 1 - 2/x2]

#### Lotka-Volterra en coordonnées canoniques

% Test avec changement de variable
function E = eulerSymplectique(f,a,b,y0,M)
h=(b-a)/M;
t = 0:1000; %time vector
dt = t(2)-t(1); %time step
tic  
Y(1,1)=y0(1);  % initialisation de p
Y(2,1)=y0(2);  % innitilisation de q
options = optimset('TolX' , 1e-6);  % option d'arret : 1e-6
for i = 1:length(t)-1  % resolution du systeme f_n+1
    r = f(Y(1,i),Y(2,i));     % appel de la fonction
    Y(1,i+1)=Y(1,i) + h*r(1); % p_n+1
    Y(:, i+1) = fsolve(@(y) Y(:, i) + h*f(Y(1,i), y(2)) - y, Y(:, i),options);
end 
hold on
plot(exp(Y(1,:)))  % plot de q
plot(exp(Y(2,:)))  % plot de p
hold off
E = Y';
end

% appel fonction avce changement de variable 
f = @(x1,x2) [0 1;-1 0]*[exp(x1)-1 ; exp(x2)-2]
E = eulerSymplectique(f,0,10,[log(2);log(5)],1000)


- On voit bien que l'on obtient exactement les même graphes avant et après **le changement de variable**. Donc notre méthode de Lotka-Volterra est correcte.