%% schema d'Euler
% ce script illustre la convergence en O(h) de la methode d'Euler.
%
close all;%ferme les figures
%
f = @(t,y) y;%definit la fonction f dans y'=f(t,y)
%f=@(t,y) y.^2;
% f(t,y)=y correspond a y'=y
%f(t,y)= y.^2 correspond a y'=y^2.
t0=0;%definit l'instant initial
tfinal=4;
%tfinal=0.9;%definit l'intervalle de temps avant singularite 
y0=1;%donnee initiale

%% trace la solutions exacte
%rentrer la solution exacte de l'edo si elle est connue
xexact=(t0:1e-2:tfinal);
yexact=y0*exp(xexact);%solution exacte exponentielle
%yexact=1./(1-xexact); %solution exacte pour y'=y^2 et y(0)=1
plot(xexact,yexact,'r');
hold on;
%
pas=[];
erreur=[];
h=0.5*(tfinal-t0);%definit le pas de depart
%
%% boucle sur les pas de temps de plus en plus petits
for k=1:10
    %initialisation
    x=(t0:h:tfinal);
    N=length(x);
    yapp=zeros(1,N);
    t = t0;
    y = y0;
    yapp(1)=y0;
    % code la methode d'Euler pour le pas h
    for j=2:N
        y = y + h*feval(f,t,y);% y(t+h) = y(t) + h*f(t,y)
        yapp(j)=y;
        t = t + h;
    end
    plot(x,yapp,'.-b');%trace la solution approchee
    legend('en rouge: solution exacte', 'en bleu: solutions approchees pour pas de temps de plus en plus petits');
    title('tapez sur une touche pour voir la convergence');

%estime l'erreur pour le pas h
    err= max(abs(yapp-exp(x)));
    %err=max(abs(yapp-1./(1-x)));
    disp('taper sur une touche pour continuer');
    pause;
    pas=[pas,h];
    erreur=[erreur,err];
    h=h/2;%divise le pas par deux.
    
end
%% graphes de l'erreur en fonction du pas.
%trace l'erreur en fonction du pas log
clf;
loglog(pas,erreur, 'b-s'); hold on;
loglog(pas, pas,'r');
legend('en bleu: erreur en norme max', 'en rouge: droite de pente 1');
title('erreur   fonction du pas en O(h), echelle log');
xlabel('pas');ylabel('erreur');
%calcule l'erreur avec le pas le plus petit.
disp('erreur finale en norme max=');disp(err);