// On résout par la méthode d'Euler explicite l'équation du pendule simple
// x"+sin(x)=0 que l'on réécrit sous la forme d'un système Y'=f(Y) où
// f(Y) = [y;sin(x)] pour Y=(x;y)

// Création d'une fonction qui retourne f(Y)

function Yout = pendule(Yin)
    Yout = [Yin(2);-sin(Yin(1))];
endfunction

// Implémentation de la méthode d'Euler explicite
// de la forme eulerexplicite(T,dt,CI)
// où T est le temps final, dt est le pas de temps
// et CI est la condition initiale

function Xsol = eulerexplicite(T,dt,CI)
    N = T/dt;
    Xsol = zeros(2,N+1);
    Xsol(:,1) = CI;
    for i=1:N
        Xsol(:,i+1)=Xsol(:,i)+dt*pendule(Xsol(:,i));
    end
endfunction

// Implémentation de la méthode de Runge-Kutta d'ordre 4
// de la forme RK4(T,dt,CI)
// où T est le temps final, dt est le pas de temps
// et CI est la condition initiale

function Xsol = RK4(T,dt,CI)
    N = T/dt;
    Xsol = zeros(2,N+1);
    Xsol(:,1) = CI;
    for i=1:N
        k1 = pendule(Xsol(:,i));
        k2 = pendule(Xsol(:,i)+(dt/2)*k1);
        k3 = pendule(Xsol(:,i)+(dt/2)*k2);
        k4 = pendule(Xsol(:,i)+dt*k3);
        Xsol(:,i+1)=Xsol(:,i)+(dt/6)*(k1+2*k2+2*k3+k4);
    end
endfunction

// Résolution numérique pour une condition initiale donnée

Tf = 100; // temps final
Delta_t = 0.05; // pas de temps
Y0 = [%pi/4;0]; // condition initiale
Solee = eulerexplicite(Tf,Delta_t,Y0); // calcul de la solution par méthode d'Euler
Solrk4 = RK4(Tf,Delta_t,Y0); // calcul de la solution par méthode RK4
t = 0:Delta_t:Tf; // discrétisation du temps

// Affichage des solutions, composante par composante, en noir Euler et en rouge RK4

figure(1)
plot(t,Solee(1,:),'black')
plot(t,Solrk4(1,:),'red')

figure(2)
plot(t,Solee(2,:),'black')
plot(t,Solrk4(2,:),'red')

// Affichage de la trajectoire dans le plan de phase
figure(3)
plot2d(Solee(1,:),Solee(2,:),rect=[-5,-2.5,5,2.5])
plot(Solrk4(1,:),Solrk4(2,:),'red')
xtitle('Trajectoire dans le plan de phase')

// Affichage de l'énergie totale du système au cours du temps
Etotee = 1-cos(Solee(1,:))+(1/2)*(Solee(2,:)).^2;
Etotrk4 = 1-cos(Solrk4(1,:))+(1/2)*(Solrk4(2,:)).^2;
figure(4)
plot(t,Etotee,'black')
plot(t,Etotrk4,'red')
xtitle('Enerige totale au cours du temps')
