%% methode des trapezes implicites
% ce script illustre la convergence en O(h^2) de la methode des trapezes
% implicites et sa stabilité par rapport à un schéma explicite plus précis
% ode45
%
close all;%ferme les figures
%% question 1
%
f = @(t,y) 3*y-4*exp(-t);%definit la fonction f dans y'=f(t,y)
%
t0=0;%definit l'instant initial
tfinal=5;%definit l'intervalle de temps 
y0=1;%donnee initiale

% trace la solutions exacte
xexact=(t0:1e-2:tfinal);
yexact=exp(-xexact);
figure(1);
plot(xexact,yexact,'red');
hold on;
%
pas=[];
erreur=[];
h=1;% pas de depart
nbiter=10; % on passe DIX fois dans la boucle
%le dernier pas de temps utilise est 1/2^9 
% la mise à jour du pas s'effectue en effet après le calcul de 
% la solution approchee

% boucle sur les pas de temps de plus en plus petits
for k=1:nbiter
    %initialisation
    x=(t0:h:tfinal);
    N=length(x);
    yapp=zeros(1,N);
    t = t0;
    y = y0;
    yapp(1)=y0;
    % code la methode des trapezes pour le pas h
    for j=2:N
        y=yapp(j-1);
        %initialisation avec Euler
        z=y+h* feval(f,t,y);
        % iterations de point fixe
        for i=1:3
            z = y + 0.5*h*(feval(f,t,y) + feval(f,t+h,z));
        end
        yapp(j)=z;
        %pas de temps suivant
        t = t + h;
    end

%estime l'erreur pour le pas h
    err=max(abs(yapp-exp(-x)));
    pas=[pas,h];
    erreur=[erreur,err];
    h=h/2;%divise le pas par deux. c'est pour celqa que nbiter=10 et pas 9
    
end
plot(x,yapp,'-k');%trace la solution approchee
legend('rouge: solution exacte', 'noir: solution approchee pour pas de temps final');
xlabel('t');
ylabel('y');

%%  Question 2
%graphes de l'erreur en fonction du pas.
%trace l'erreur en fonction du pas log
figure(2);% une nouvelle figure
loglog(pas,erreur, 'blue -s'); hold on;
loglog(pas, pas.^2,'red -s');
legend('en bleu: erreur en norme max', 'en rouge: droite de pente 2');
title('erreur  en O(h^2), echelle log');
xlabel('pas');ylabel('erreur');
%calcule  et affiche l'erreur avec le pas le plus petit.
disp('erreur finale en norme max=');disp(err);
% Le schéma est d'ordre deux car en loglog  les erreurs sont
% sur une droite de pente 2
%% Question 3
%solution standard matlab
figure(3);
plot(xexact,yexact,'r')
hold on;
ode45(@(t,y) 3*y-4*exp(-t),[0,4.2],1) % tfinal =4.1 pour voir ce qui se passe apres
legend(' solution exacte', 'solution matlab ode45');
xlabel('t');
ylabel('y');
hold off;

%%% le fait que Matlab ode45 ne marche pas bien vient du fait que la 
%%%% solution generale de l'EDO est C*exp(3*t) + exp(-t). Comme les calculs
%%%% numériques sont effectués avec une précision finie, les valeurs
%%%% calculées correspondent à une solution  dans laquelle
%%% la constante C=eps n'etant pas parfaitement nulle et est de l'ordre de la precision machine
%%% la composante eps*exp(3*t) augmente exponentiellement au cours du temps
%%% et finit par dominer la solution exacte. Si on augmente trop tfinal la
%%% composante C*exp(3*t) devient trop grande et masque complètement la
%%% solution exacte. 
%%% La méthode des trapèzes implicite va aussi s'écarter peu à peu de la solution
%%% exacte bien qu'elle semble un peu meilleure que ode45 ici
%% tracé de solutions voisines de la solution qui s'en écartent progressivement.
figure(4);
xexact=(t0:1e-2:tfinal);
plot(xexact, exp(-xexact),'r'); hold on
for C = [-3e-7  -2e-7  -1e-7 -0.5e-7 0.5e-7  1e-7  2e-7 3e-7 ]
    plot(xexact, C*exp(3*xexact)+ exp(-xexact))
    hold on
end
xlabel('t')
ylabel('y')
title('solutions Cexp(3t) + exp(-t) pour -3e-7< C < 3e-7')