// On résout par la méthode d'Euler implicite l'équation différentielle suivante
// x'(t) = 1-cos(x(t)) où l'on connait les solutions exactes et ainsi on pourra
// estimer l'erreur relative commise par la méthode d'Euler implicite et retrouver
// que cette dernière est de l'ordre du pas de temps

// on commence par implémenter une méthode de Newton pour résoudre
// z-af(z)-b = 0 où f(z)=1-cos(z)

function z = newton(a,b)
    z = b;
    for k=1:10
        z = z-(z-a*(1-cos(z))-b)/(1-a*sin(z));
    end
endfunction

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

function Xsol = eulerimplicite(T,dt,CI)
    N = T/dt;
    Xsol = zeros(1,N+1);
    Xsol(1,1) = CI;
    for i=1:N
        Xsol(1,i+1) = newton(dt,Xsol(1,i));
    end
endfunction

// Résolution numérique pour une condition initiale donnée
// et comparaison avec la solution exacte

Tf = 10; // temps final
Delta_t = 0.1; // pas de temps
x0 = %pi; // condition initiale
Solapp = eulerimplicite(Tf,Delta_t,x0); // calcul de la solution par méthode d'Euler
t = 0:Delta_t:Tf; // discrétisation du temps 
tt = Delta_t:Delta_t:Tf;
Solex = [%pi 2*atan(-1 ./tt)+2*%pi];

figure(1)
plot(t,Solapp,'black')
plot(t,2*%pi*ones(1,Tf/Delta_t+1),'red')
plot(t,Solex,'blue')

// Calcul de l'erreur de la méthode d'Euler implicite en fonction du pas de temps

T = 20;
ci = %pi;
dt = 0.005:0.005:0.1;
m = length(dt);
Erreur = zeros(1,m);
for k=1:m
    Dt = dt(k);
    Tsol = Dt:Dt:T;
    Solexact = [%pi 2*atan(-1 ./Tsol)+2*%pi];
    Erreur(k) = max(abs(eulerimplicite(T,Dt,ci)-Solexact))
end

// Permet de calculer l'équation de la  "meilleure" (au sens des moindres carrés) droite 
// passant par un ensemble de points donnés

function cf = polyfit(x,y,n)
A = ones(length(x),n+1)
for i=1:n
    A(:,i+1) = x(:).^i
end
cf = lsq(A,y(:))
endfunction

Y = polyfit(dt,Erreur,1);


// Visualisation de l'erreur et on super-impose la droite interpolant cette erreur
// ce qui permet de confirmer que l'erreur est linéaire en le pas de temps


figure(2)
plot(dt,Erreur,'d')
plot(dt,Y(1)+dt*Y(2),'red')
xtitle('Erreur de la méthode Euler implicite')













