lab5

clc;
close all;
format long e;

%Question 1
%a)

%  y'(t) = 1 + y(t) / t
%  y'(t) - y(t) / t = 1
%
%  Facteur intégrant : u(t) = exp(integral(-dt / t))
%                           = exp(-ln(t))
%                           = t ^ -1 = 1 / t
% y(t) = t * (integral(u(t) * dt) + c)
%      = t (ln(t) + c)
%      = t * ln(t) + t*c
% y(1) = 1 * ln(1) + 1 * c
% 1 = 1 * c
% 1 = c
%
% Réponse : y(t) = t * ln(t) + t
%                = t (ln(t) + 1)

%b)
t0 = 1;
y0 = 1;
h = 0.1;
nbPas = 50;

[t,yQ1] = ptmilieu('prob1', t0, y0, h, nbPas);

yt = t .* (log(t) +1);

figure(1);
plot(t, yQ1, 'b--+', t, yt, 'm-o');
title('Méthode de point milieu vs solution analytique');
xlabel('t');
ylabel('y(t)');
legend('Point Milieu', 'Solution analytique');

%c)
yPtMil = [1:51];
for i = 1:51
    yPtMil(i) = yQ1(i);
end

%Erreur relative max:
errRelMax = max(abs(yPtMil-yt)./yt);
fprintf('Erreur relarive max: %.10e \n',errRelMax)


eh = 1:1:5;

fprintf('\n\n----------------------------------------------------------\n')
fprintf('h \t\t\t\t e(h) \t\t\t\t\t e(h)/e(h/2)\n')
fprintf('----------------------------------------------------------\n')
for i = 1:5

    h2 = 1/2^(i+1);
    nbPas2 = 5 ./h2;
    [t,yQ1c] = ptmilieu('prob1', t0, y0, h2, nbPas2);

    %figure(i+1);
    %plot(t, yQ1c, 'b--+');
    %title('Méthode de point milieu vs solution analytique');
    %xlabel('t');
    %ylabel('y(t)');
    %legend('Point Milieu');
    %fprintf('estimé : %f \n',yQ1c(nbPas2+1))
    eh(i) = yQ1c(nbPas2+1)-yt(nbPas+1);

    if(i > 1)
        fprintf('%f \t\t %.10e \t\t %.10e\n',h2,eh(i),eh(i-1)/eh(i))
    else
        fprintf('%f \t\t %.10e \t\t ----\n',h2,eh(i))
    end

end
fprintf('----------------------------------------------------------\n')

%L'ordre de convergence est donc de 2 puisqu'on a :
% 2^ordre = e(h)/e(h/2) approximativement = 4

%Question 2
%a)

%y(t) = 5 e^(-7 t)


%b)
erreurAbs1 = zeros((2/0.5)+1);
erreurAbs2 = zeros((2/0.25)+1);
erreurAbs3 = zeros((2/0.125)+1);



for i = 1:3

    h3 = 1/(2^i);
    nbPas3 = 2 / h3;
    [t3, yEuler] = euler('prob2', 0, 5, h3, nbPas3);

    yAnal = 5 * exp(-7 * t3);

    yEuler2 = [1:nbPas3+1];
    for j = 1:nbPas3+1
        yEuler2(j) = yEuler(j);
    end

    switch i
        case 1
            erreurAbs1 = abs(yAnal-yEuler2);
        case 2
            erreurAbs2 = abs(yAnal-yEuler2);
        case 3
            erreurAbs3 = abs(yAnal-yEuler2);
    end



    figure(i + 1);
    plot(t3, yEuler, 'b--+', t3, yAnal, 'm-o');
    title(sprintf('Méthode euler h=%f vs solution analytique', h3));
    xlabel('t');
    ylabel('y(t)');
    legend('Euler', 'Solution analytique');

end

t1b = 0:0.5:2;
t2b = 0:0.25:2;
t3b = 0:0.125:2;

figure(5);
plot(t1b, erreurAbs1, 'b--+', t2b, erreurAbs2, 'm--*', t3b, erreurAbs3, 'k--x');
title('Erreurs absolues dans le calcul par Euler selon le pas h');
xlabel('t');
ylabel('Erreur absolue');
legend('h=0.5', 'h=0.25', 'h=0.125');

%c)

%Il faut que le pas soit le plus près possible de 0 tout en restant loin de
%l'erreur machine.

a = 1;
h_exp = 0.25; %on sait de par le point b) qu'avec 0.125 la solution d'euler tend vers 0
%on sait également qu'avec 0.25 euler ne tend plus vers 0 (de facon stritement décroissante)
fprintf('\n----------------------------------------------------------\n')
fprintf('Vérification de la décroissance selon les h\n');
fprintf('----------------------------------------------------------\n')
for h = 0.125:0.01:0.25
    nbpas = 2/h;
    [t, y] = euler('prob2', t0, y0, h, nbpas);
    yAnal = 5 * exp(-7 * t);

    %on vérifie si trois points consécutifs décroissent
    decroissant = 1;
    if(y(1) - y(2) < 0)
        decroissant = 0;
    end

    if(y(2) - y(3) < 0)
        decroissant = 0;
    end

    %On peut également vérifier si euler converge vers la solution analytique
    %en regardant si l'erreur relative rétressie

    retressie = 1;
    for j = 1:nbpas

        errRel = (y(j) - yAnal(j))/yAnal(j);
        errRelNext = (y(j+1) - yAnal(j+1))/yAnal(j+1);

        if(errRel < errRelNext)
            retressie = 0;
        end

    end

    if(retressie == 1)
       fprintf('h = %f\tconverge\n',h);
    else
       fprintf('h = %f\tne converge pas\n',h);
    end

    if(decroissant == 1)
       fprintf('\t\t\t\tstritement décroissant\n\n',h);
    else
       fprintf('\t\t\t\toscille\n\n',h);
    end
end

%d(
%On semble avoir euler explicite alors on doit référer à la page 46 (49)
%des notes de cours comlémentaires

%Par la méthode donnée à cette page, nous avons pu démontrer que la
%méthode d'Euler est stable pour un h se situant entre 0 et 2/7