%% convergence des schemas
% ce script illustre la convergence en O(h) de la methode d'Euler
%la convergence en O(h^2) de la methode de Runge du point milieu
%et la convergence en O(h^3) de la methode de Heun
%
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=7;%definit l'intervalle de temps 
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,'black');
hold on;

pas=[];
erreur_euler=[];
erreur_heun=[];
erreur_milieu=[];
erreur_rk4=[];
h=1;%definit le pas de d�part
%
%% boucle sur les pas de temps de plus en plus petit
for k=1:5
    %initialisations
    x=(t0:h:tfinal);
    N=length(x);
    yeuler=zeros(1,N);
    t = t0;
    y = y0;
    yeuler(1)=y0;
    %methode d'Euler
    for j=2:N
        y = y + h*feval(f,t,y);
        yeuler(j)=y;
        t = t + h;
    end
    %methode du point milieu
    ymilieu=zeros(1,N);
    t=t0;
    y=y0;
    ymilieu(1)=y0;
    for j=2:N
        u = y + (h/2)*feval(f,t,y);
        y=y+ h*feval(f,t+h/2,u);
        ymilieu(j)=y;
        t = t + h;
    end
    %methode du trapeze explicite RK2
    %for j=2:N
    %    u = y + h*feval(f,t,y);
    %    y=y+ h*1/2*(feval(f,t,y)+feval(f,t+h,u));
    %    yapp(j)=y;
    %   t = t + h;
    %end
    %methode de Heun
    yheun=zeros(1,N);
    t=t0;
    y=y0;
    yheun(1)=y0;
    for j=2:N
        u2 = y + h/3*feval(f,t,y);
        u3= y + 2/3*h*feval(f,t+h/3,u2);
        y=y+ h*(1/4*feval(f,t,y)+3/4*feval(f,t+2*h/3,u3));
        yheun(j)=y;
        t = t + h;
    end
    %methode de Runge Kutta 4
    yrk4=zeros(1,N);
    t=t0;
    y=y0;
    yrk4(1)=y0;
    for j=2:N
        u2 = y + h/2*feval(f,t,y);
        u3= y + h/2*feval(f,t+h/2,u2);
        u4= y + h*feval(f,t+h/2,u3);
        y=y+ h*( 1/6*feval(f,t,y)+ 2/6*feval(f,t+h/2, u2)+ 2/6*feval(f,t+h/2,u3) + 1/6*feval(f,t+h,u4) );
        yrk4(j)=y;
        t = t + h;
    end
    % trace des courbes pour le pas le plus grossier
    if k==1 
        plot(x,yeuler,'red -s',x,ymilieu,'green-^', x,yheun,'blue -o',x,yrk4,'m -x');%trace les solution approchee
        legend('noir: solution exacte', 'rouge: schema Euler','vert: schema  Runge point milieu', 'bleu: schema Heun','magenta: Runge Kutta 4','Location','NW');
        title('Les differents schemas explicites a un pas');
        xlabel('t');
        ylabel('y(t)');
    end;

%estime l'erreur pour le pas h
    err_euler= max(abs(yeuler-exp(x)));
    err_milieu=max(abs(ymilieu-exp(x)));
    err_heun= max(abs(yheun-exp(x)));
    err_rk4= max(abs(yrk4-exp(x)));
    %err=max(abs(yapp-1./(1-x)));
    
    pas=[pas,h];
    erreur_euler=[erreur_euler,err_euler];
    erreur_milieu=[erreur_milieu,err_milieu];
    erreur_heun=[erreur_heun,err_heun];
    erreur_rk4=[erreur_rk4,err_rk4];
    %%%%%%%%%%
    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 
figure; % cree une nouvelle figure
loglog(pas,erreur_euler, 'red -s'); hold on;
loglog(pas,erreur_milieu, 'green-^');
loglog(pas,erreur_heun, 'blue -o');
loglog(pas,erreur_rk4, 'm -x');
% %trace des droites de pente 1,2,3,4 en log log
% loglog(pas, pas.^4,'m:');
% loglog(pas, pas.^3,'b:');
% loglog(pas, pas.^2,'g:');
% loglog(pas, pas,'r:');
legend('rouge: erreur en norme max schema Euler',...
     'vert: schema Runge point milieu','bleu: schema de Heun',...
     'magenta: RK4','Location', 'SE','AutoUpdate','off');
 title('erreur en  fonction du pas, echelle log log');
xlabel('pas');ylabel('erreur');
%%
%pente exacte des courbes
%trace des triangles de pente
%affichage des pentes exactes
line([pas(5),pas(4)],[erreur_euler(5), erreur_euler(5)],'Color','r');
line([pas(4),pas(4)],[erreur_euler(5), erreur_euler(4)],'Color','r');
str={['pente=' num2str( (log10(erreur_euler(5))-log10(erreur_euler(4)))/(log10(pas(5))-log10(pas(4))),3 )]};
text(pas(4), (erreur_euler(5)+erreur_euler(4))/2,str)
%
line([pas(5),pas(4)],[erreur_milieu(5), erreur_milieu(5)],'Color','g');
line([pas(4),pas(4)],[erreur_milieu(5), erreur_milieu(4)],'Color','g');
str={['pente='...
    num2str(...
    (log10(erreur_milieu(5))-log10(erreur_milieu(4)))...
    /(log10(pas(5))-log10(pas(4))),3 )]};
text(pas(4), (erreur_milieu(5)+erreur_milieu(4))/2,str)
%
line([pas(5),pas(4)],[erreur_heun(5), erreur_heun(5)],'Color','b');
line([pas(4),pas(4)],[erreur_heun(5), erreur_heun(4)],'Color','b');
str={['pente=' num2str( (log10(erreur_heun(5))-log10(erreur_heun(4)))/(log10(pas(5))-log10(pas(4))),3 )]};
text(pas(4), (erreur_heun(5)+erreur_heun(4))/2,str);
%
line([pas(5),pas(4)],[erreur_rk4(5), erreur_rk4(5)],'Color','m');
line([pas(4),pas(4)],[erreur_rk4(5), erreur_rk4(4)],'Color','m');
str={['pente=' ...
    num2str(...
    (log10(erreur_rk4(5))-log10(erreur_rk4(4)))/(log10(pas(5))-log10(pas(4))),3 )]};
text(pas(4), (erreur_rk4(5)+erreur_rk4(4))/2,str)
%
%%calcule l'erreur avec le pas le plus petit.
%
disp('erreurs finales en norme max=');err_euler,err_milieu,err_heun,err_rk4
disp('le pas de temps est');h