function LevenbergMarquardtMethod

%p est un vecteur de dimension 2x1 contenant les parametres

p(1,1) = 0.4;
p(2,1) = 0.2;

%F = y - (p(1).*t)./(p(2)+t); -> forme de l'equation a minimiser

%dF1 = -t./(p(2)+t);          -> forme de la premiere derivee en p(1)
%dF2 = p(1)*x./((p(2)+t).^2); -> forme de la seconde derivee en p(2)

%donnees experimentales
t = [0.038 ; 0.194 ; 0.425 ; 0.626 ; 1.253 ; 2.500 ; 3.740];
D = [0.050 ; 0.127 ; 0.094 ; 0.2122 ; 0.2729 ; 0.2665 ; 0.3317];

%initialisation de lambda
lambda0 = 0.01;
rmsea = 15;
rmse = 10;
inc = 0;

%for inc = 1:10
while abs(rmsea-rmse(end))>1e-5 && inc<500 && lambda0 < 1e7
    
    inc = inc+1;
    
    %calcul du Jacobien
    J(1:length(t),1) = -t./(p(2)+t);
    J(1:length(t),2) = p(1)*t./((p(2)+t).^2);
    
    %calcul des residus
    F = D-(p(1).*t)./(p(2)+t);
    
    %calul de la rmse lors de la premiere iteration
    if inc==1
        rmse = (sum(F.^2)/length(F))^0.5;
    else
        %mise a jour de lambda
        lambda = lambda0;
        %generation de la matrice unite
        tailleJJ = length(J.'*J);
        %calul des delta en fonction de lambda
        Delta1 = -inv((J.'*J)+lambda*eye(tailleJJ))*(J.')*F;
        lambda=lambda0/10;
        Delta2 = -inv((J.'*J)+lambda*eye(tailleJJ))*(J.')*F;
        p1 = p;
        p2 = p;
        %mise a jour des parametres dans un nouveau vecteur
        p1 = p1+Delta1;
        %calcul des nouveaux residus et de la nouvelle
        %rmse pour lambda=lambda0
        F = D-(p1(1).*t)./(p1(2)+t);
        rmse1 = (sum(F.^2)/length(F))^0.5;
        %calcul des nouveaux residus et de la nouvelle
        %rmse pour lambda=lambda0/10
        p2 = p2+Delta2;
        F = D-(p2(1).*t)./(p2(2)+t);
        rmse2 = (sum(F.^2)/length(F))^0.5;
        %condition sur la mise a jours des parametres
        %voir texte pour plus de details
        if rmse1>rmse && rmse2>rmse
            lambda0 = lambda0*10;
        elseif rmse2<rmse1
            lambda0 = lambda0/10;
            p = p2;
            rmsea = rmse;
            rmse = rmse2;
        elseif rmse2>=rmse1
            p = p1;
            rmsea = rmse;
            rmse = rmse1;
        end
        
    end
    
end

%calcul de la rmse
rmse = (sum(F.^2)/length(F))^0.5;

%calcul de l'ecart type (standard deviation)
stdv = (sum(F.^2)/(length(F)-length(p)))^0.5;

%evalutation des donnees ajustees
ta = 0:0.01:4;
Da = (p(1).*ta)./(p(2)+ta);

%affichage des resultats
fprintf('\np1 = %f\ntp2 = %f\n\n',p);
fprintf('rmse = %f\n',rmse);
fprintf('stdv = %f\n\n',stdv);

%affichage du graphique
figure

h = plot(t, D, 'or', ta, Da, 'b');

legend(h, {'Donnees experimentales' 'Meilleur ajustement'}, 'location', 'southeast')

title('Ajustement avec la methode de Levenberg-Marquardt');
xlabel('t');
ylabel('Donnes a ajuster');