// ###################################################
// ######### Laurent Dietrich ########################
// ## http://math.univ-toulouse.fr/~ldietric/ ########
// ###################################################
// Ceci est un exemple de code plus ou moins complet
// sur un texte d'équations différentielles de type
// agrégation (option B). Toute remarque sera la bien-
// venue, cf. ma page perso pour me contacter.
// ###################################################



clear; // efface toute la mémoire
xdel(winsid()); // efface les figures
clc; // nettoie l'invite de commande

// paramètres du modèle
a = 2;
b = 3;
c = 1;
d = 4;


// ###################  On va définir ici toutes nos fonctions ################
// Fonction F .................................................................
deff('v = F(u)', 'v(1)=(a-b.*u(2)).*u(1);v(2)=(-c+d.*u(1)).*u(2);');
// ............................................................................

// Fonction H .................................................................
deff('t = H(x,y)', 't = d*x-c*log(x)+b*y-a*log(y)');
// ............................................................................


// Fonction Euler exp ........................................................
function [u]=eulerexp(u0,dt,T) // la fonction eulerexp demande une position initiale, un pas de discrétisation, et un temps d'arrêt
   
    u = [u0];
    N = floor(T/dt); // le nombre de points pour arriver jusqu'à (environ) T : je préfère choisir dt et en déduire N, c'est plus mathématique (une méthode converge quand dt -> 0 donc c'est bien de pouvoir le contrôler)
    
    for i=1:N
         uold = u(:,$); // on récupère la dernière colonne
         u = [u, uold + dt*F(uold)]; // on concatène en colonne des couples de ligne
    end
endfunction
// ...........................................................................

// Fonction RK4 ........................................................
function [u]=rk4(u0,dt,T)
   
    u = [u0];
    N = floor(T/dt);
    
    for i=1:N
         uold = u(:,$);
         k1 = uold;
         k2 = uold + dt/2*F(k1);
         k3 = uold + dt/2*F(k2);
         k4 = uold + dt*F(k3);
         u = [u, uold + dt/6*(F(k1)+2*F(k2)+2*F(k3)+F(k4))];
    end
endfunction
// ...........................................................................

// Fonction Heun ........................................................
function [u]=heun(u0,dt,T)
   
    u = [u0];
    N = floor(T/dt);
    
    for i=1:N
         uold = u(:,$);
         uinter = uold + dt*F(uold);
         u = [u, uold + dt/2*(F(uold)+F(uinter))];
    end
endfunction
// ...........................................................................
// ###########################################################################




// ########## Ici commence le script avec les éxécutions ####################

// .... Dessin du champ de vecteurs ..........................................
N = 20; // taille de la grille pour le dessin du c.v.
x = linspace(0,2*c/d,N);
y = linspace(0,2*a/b,N);
[X,Y] = meshgrid(x,y);

for i = 1:N
    for j = 1:N
        cx(i,j) = (a-b*Y(j,i)).*X(j,i); // on met des j,i car, cf. le [X,Y], la valeur de x évolue sur les COLONNES de X et celle de y sur les LIGNES de Y !
        cy(i,j) = (-c+d*X(j,i)).*Y(j,i);
    end
end

champ1(x,y,cx,cy);
cv = gce();
cv.thickness = 0.6;
title('Champ de vecteurs'); // mettre un titre au graphe
xlabel('x'); // donner un nom à l'axe des 
ylabel('y');
// ...........................................................................


halt("Appuyez sur entrée"); // on demande une pause jusqu'à appuyer entrée dans la commande
u0 = [0.4;0.9]; // Donnée initiale commune à tout ce qui suit

// Calcul via Euler exp ...............................................................
T = 50;
dt = 0.01;
t = 0:dt:T;

u = eulerexp(u0,dt,T);
x = u(1,:);
y = u(2,:);
figure; // ouvre une nouvelle figure
subplot(311); // la déclare comme un tableau 3x1, et on tracera sur la première (en haut)
plot(t,x,t,y);
title('Évolution des densités - Euler exp. T=50, h=0.01');
xlabel('t');
legend('x(t)','y(t)');

for i = 1:length(t) // on retient l'énergie numérique
    EE(i) = H(x(i),y(i));
end
// ...........................................................................

// RK4 ...............................................................
T = 50;
dt = 0.1;
tr = 0:dt:T;
ur = rk4(u0,dt,T);
xr = ur(1,:);
yr = ur(2,:);
subplot(312); // on trace sur la deuxième figure (milieu)
plot(tr,xr,tr,yr);
title('Évolution des densités - RK4. T=50, h=0.1');
xlabel('t');
legend('x(t)','y(t)');

for i = 1:length(tr)
    ER(i) = H(xr(i),yr(i));
end
// ...........................................................................

// Heun ...............................................................
T = 50;
dt = 0.1;
th = 0:dt:T;
uh = heun(u0,dt,T);
xh = uh(1,:);
yh = uh(2,:);
subplot(313); // on tracera sur la 3ème figure en bas
plot(th,xh,th,yh);
title('Évolution des densités - Heun. T=50, h=0.1');
xlabel('t');
legend('x(t)','y(t)');

for i = 1:length(tr)
    Eh(i) = H(xh(i),yh(i));
end
// ...........................................................................


halt("Appuyez sur entrée");

// Plan de phase .............................................................
figure; // nouvelle figure

subplot(311);
plot(x,y,'b',u0(1),u0(2),'rx');
xlabel('x(t)');
ylabel('y(t)');
title('Évolution dans le plan de phase - Euler exp.');

subplot(312);
plot(xr,yr,'b',u0(1),u0(2),'rx');
xlabel('x(t)');
ylabel('y(t)');
title('Évolution dans le plan de phase - RK4');

subplot(313);
plot(xh,yh,'b',u0(1),u0(2),'rx');
xlabel('x(t)');
ylabel('y(t)');
title('Évolution dans le plan de phase - Heun');

// ...........................................................................

halt("Enteeeer ?");

// Énergie ...................................................................
figure;
subplot(311);
plot(t,EE);
title('Hamiltonien H - Euler exp.');
xlabel('t');
ylabel('H(u(t))');

subplot(312);
plot(tr,ER);
title('Hamiltonien H - RK4');
xlabel('t');
ylabel('H(u(t))');

subplot(313);
plot(th,Eh);
title('Hamiltonien H - Heun');
xlabel('t');
ylabel('H(u(t))');


// ################### Ici on code et on illustre Euler implicite #############
// ... Fonctions g et Dg pour le Newton ......................................
function t = g(x,un,h)
    t = x - un - h*F(x); // résoudre g(x) = 0 correspond à résoudre l'équation x = un + h*F(x). Si F est lipschitzienne et h assez petit, g(x)-x sera contractante.
endfunction

function M = Dg(x,h) // la différentielle de g d'avant, calculé explicitement puisqu'on peut le faire ici ! (Sinon, le faire numériquement...)
    M = eye(2,2) - h*[a-b*x(2),-b*x(1);d*x(2),-c+d*x(1)]; // qui sera bien sûr inversible si h est assez petit
endfunction
// ...........................................................................

// Fonction Newton ........................................................
function [u] = newton(un,h,N)
    u = un; // POINT DE DÉPART de Newton (ie le x_0): il doit être proche de la solution, on lui donnera donc le un, ie le point précédent dans le calcul de la sol numérique, en espérant que c'est suffisant. Sinon, on peut aussi partir de la sol. d'Euler explicite, c'est mieux ! Ici ça suffira.
    for i=1:N
        u = u - Dg(u,h)\g(u,un,h);        
    end
endfunction
// ...........................................................................

// Fonction Euler imp ........................................................
function [u]=eulerimp(u0,dt,T) // Tout comme avant sauf qu'il faut sortir un+1 via un Newton qu'on vient de faire :)
   
    u = [u0];
    N = floor(T/dt);
    
    for i=1:N
         uold = u(:,$);
         u = [u, newton(uold,dt,5)];
    end
endfunction
// ...........................................................................

halt('Entrée pour passer à l`implicite');

// Affichage Euler imp ...............................................................
T = 50;
dt = 0.01;
t = 0:dt:T;

u = eulerimp(u0,dt,T);
x = u(1,:);
y = u(2,:);
figure;
plot(t,x,t,y);
title('Évolution des densités - Euler imp. T=50, h=0.01');
xlabel('t');
legend('x(t)','y(t)');

for i = 1:length(t)
    EI(i) = H(x(i),y(i));
end

halt('Enter');

figure;

plot(x,y,'b',u0(1),u0(2),'rx');
xlabel('x(t)');
ylabel('y(t)');
title('Évolution dans le plan de phase - Euler imp.');

halt('Enter');

figure;
plot(t,EI);
title('Hamiltonien H - Euler imp.');
xlabel('t');
ylabel('H(u(t))');


// ...........................................................................


