207 lines
7.2 KiB
Matlab
Executable file
207 lines
7.2 KiB
Matlab
Executable file
%-------------------------------------------------------------------------%
|
||
% 1SN - TP Optimisation %
|
||
% INP Toulouse - ENSEEIHT %
|
||
% Novembre 2020 %
|
||
% %
|
||
% Ce fichier contient le programme principal permettant l'estimation %
|
||
% des parametres de la fonction de desintegration radioactive du %
|
||
% carbone 14 par une approche des moindres carres. % %
|
||
% Modele : A(t)= A0*exp(-lambda*t) %
|
||
% Les algorithmes utilises pour la minimisation sont : %
|
||
% - l'algorithme de Gauss-Newton %
|
||
% - l'algorithme de Newton. %
|
||
%-------------------------------------------------------------------------%
|
||
|
||
clear
|
||
close all
|
||
clc
|
||
format shortG
|
||
taille_ecran = get(0,'ScreenSize');
|
||
L = taille_ecran(3);
|
||
H = taille_ecran(4);
|
||
|
||
% Initialisation
|
||
|
||
% Donnees
|
||
if exist('C14_results.txt','file')
|
||
delete('C14_results.txt');
|
||
end
|
||
diary C14_results.txt
|
||
Ti = [ 500; 1000; 2000; 3000; 4000; 5000; 6300];
|
||
Ai = [14.5; 13.5; 12.0; 10.8; 9.9; 8.9; 8.0];
|
||
Donnees = [Ti, Ai];
|
||
|
||
% Estimation a priori des parametres du modele : beta0 = [A0, lambda]
|
||
beta0 = [10; 0.0001]; % Newton, Gauus-Newton, fminunc et leastsq converge
|
||
% beta0 = [15; 0.001]; % Newton, Gauss-Newton, fminunc et leastsq divergent
|
||
% beta0 = [15; 0.0005]; % Newton diverge, Gauss-Newton, fminunc et leastsq convergent
|
||
% beta0 = [10; 0.0005]; % Gauss-Newton converge
|
||
|
||
%% Calcul et affichage du modele initial ----------------------------------
|
||
xmin = 9; xmax = 20;
|
||
x = linspace(xmin, xmax, 100);
|
||
|
||
ymin = -0.0001; ymax = 0.0005;
|
||
y = linspace(ymin, ymax, 100);
|
||
|
||
[A0_plot, lambda_plot] = meshgrid(x, y);
|
||
[m , n] = size(A0_plot);
|
||
f_plot = zeros(m, n);
|
||
|
||
for i = 1:m
|
||
for j = 1:n
|
||
res_plot = residu_C14([A0_plot(i,j); lambda_plot(i,j)], Donnees);
|
||
f_plot(i,j) = res_plot.'*res_plot /2;
|
||
end
|
||
end
|
||
|
||
figure('Position',[0.1*L,0.1*H,0.8*L,0.8*H]);
|
||
subplot(1,3,1) % Fonction f en 3D
|
||
mesh(A0_plot, lambda_plot, f_plot)
|
||
axis('square')
|
||
title('Representation de la fonction f des moindres carres')
|
||
xlabel('A_0');
|
||
ylabel('\lambda')
|
||
zlabel('f(A_0,\lambda)')
|
||
|
||
T = linspace(0,6500,100);
|
||
A = beta0(1)*exp(-beta0(2)*T);
|
||
txt_legend{1} = 'donnees';
|
||
txt_legend{2} = 'depart';
|
||
|
||
subplot(2, 3, 2) % Donnees + modele initial pour Gauss-Newton
|
||
title('Desintegration radioactive du carbone 14 (Gauss-Newton)')
|
||
axis([0 max(T) 0 18])
|
||
xlabel('dur<EFBFBD>e T')
|
||
ylabel('radioactivit<EFBFBD> A')
|
||
hold on
|
||
grid on
|
||
plot(Ti, Ai, 'ok')
|
||
plot(T, A)
|
||
legend(txt_legend{1:2}, 'Location', 'SouthWest')
|
||
|
||
subplot(2, 3, 5) % Courbes de niveaux de f pour Gauss-Newton
|
||
title('Recherche des parametres (Gauss-Newton)')
|
||
xlabel('A_0')
|
||
ylabel('\lambda')
|
||
hold on
|
||
contour(A0_plot, lambda_plot, f_plot, 100);
|
||
plot(beta0(1),beta0(2),'ok')
|
||
text(beta0(1),beta0(2),' depart \beta^{(0)}')
|
||
|
||
subplot(2, 3, 3) % Donnees + modele initial pour Newton
|
||
title('Desintegration radioactive du carbone 14 (Newton)')
|
||
axis([0 max(T) 0 18])
|
||
xlabel('dur<EFBFBD>e T')
|
||
ylabel('radioactivit<EFBFBD> A')
|
||
hold on
|
||
grid on
|
||
plot(Ti, Ai, 'ok')
|
||
plot(T, A)
|
||
legend(txt_legend{1:2}, 'Location', 'SouthWest')
|
||
|
||
subplot(2, 3, 6) % Courbes de niveaux de f pour Newton
|
||
title('Recherche des parametres (Newton)')
|
||
xlabel('A_0')
|
||
ylabel('\lambda')
|
||
hold on
|
||
contour(A0_plot, lambda_plot, f_plot, 100);
|
||
plot(beta0(1),beta0(2),'ok')
|
||
text(beta0(1),beta0(2),' depart \beta^{(0)}')
|
||
|
||
pause(0.5)
|
||
|
||
%% Algorithmes
|
||
% Choix du nombre d'iteration maximal
|
||
|
||
nb_iterations_max = 8;
|
||
txt_legend = cell(1,nb_iterations_max+2);
|
||
txt_legend{1} = 'donnees';
|
||
txt_legend{2} = 'depart';
|
||
|
||
|
||
%% Gauss-Newton
|
||
% -------------
|
||
% Initialisation de l'affichage
|
||
disp('Algorithme de Gauss-Newton')
|
||
disp('--------------------------------------------------------------------------------------------')
|
||
disp(' nb_iter A0 lambda ||f''(beta)|| f(beta) ||delta|| exitflag ')
|
||
disp('--------------------------------------------------------------------------------------------')
|
||
|
||
% Calcul et affichage des valeurs initiales
|
||
res_beta = residu_C14(beta0, Donnees);
|
||
f_beta = 0.5*(res_beta.')*res_beta;
|
||
J_res_beta = J_residu_C14(beta0, Donnees);
|
||
norm_grad_f_beta = norm((J_res_beta.')*res_beta);
|
||
disp([0 beta0(1) beta0(2) norm_grad_f_beta f_beta]);
|
||
|
||
options = [sqrt(eps) 1.e-12 0];
|
||
|
||
for i = 1:nb_iterations_max
|
||
txt_legend{i+2} = ['iteration ' num2str(i)];
|
||
options(3) = i;
|
||
[beta, norm_grad_f_beta, f_beta, norm_delta, k, exitflag] = ...
|
||
Algo_Gauss_Newton(@(beta) residu_C14(beta, Donnees), ...
|
||
@(beta) J_residu_C14(beta, Donnees), ...
|
||
beta0, options);
|
||
disp([k beta(1) beta(2) norm_grad_f_beta f_beta norm_delta exitflag])
|
||
% Visualisation
|
||
A = beta(1)*exp(-beta(2)*T);
|
||
subplot(2, 3, 2)
|
||
plot(T,A)
|
||
legend(txt_legend{1:i+2})
|
||
% eval(['print -depsc fig_GN_courbe' int2str(i) '_C14'])
|
||
subplot(2, 3, 5)
|
||
plot(beta(1),beta(2),'ok')
|
||
text(beta(1),beta(2),[' \beta^{(' num2str(i) ')}'])
|
||
|
||
pause(0.5)
|
||
end
|
||
|
||
% Affichage des itérés de beta et sauvegarde des graphique
|
||
disp('--------------------------------------------------------------------------------------------')
|
||
|
||
%% Newton
|
||
% -------
|
||
% Initialisation de l'affichage
|
||
disp('Algorithme de Newton')
|
||
disp('--------------------------------------------------------------------------------------------')
|
||
disp(' nb_iter A0 lambda ||f''(beta)|| f(beta) ||delta|| exitflag ')
|
||
disp('--------------------------------------------------------------------------------------------')
|
||
|
||
% Calcul et affichage des valeurs initiales
|
||
res_beta = residu_C14(beta0, Donnees);
|
||
f_beta = 0.5*(res_beta.')*res_beta;
|
||
J_res_beta = J_residu_C14(beta0, Donnees);
|
||
norm_grad_f_beta = norm((J_res_beta.')*res_beta);
|
||
disp([0 beta0(1) beta0(2) norm_grad_f_beta f_beta]);
|
||
|
||
options = [sqrt(eps) 1.e-12 0];
|
||
|
||
for i = 1:nb_iterations_max
|
||
options(3) = i;
|
||
% Algorithme de Newton
|
||
[beta, norm_grad_f_beta, f_beta, norm_delta, k, exitflag] = ...
|
||
Algo_Newton(@(beta) Hess_f_C14(beta, Donnees, ...
|
||
@(beta) residu_C14(beta, Donnees), ...
|
||
@(beta) J_residu_C14(beta, Donnees)), ...
|
||
beta0, options);
|
||
% Affichage des valeurs
|
||
disp([k beta(1) beta(2) norm_grad_f_beta f_beta norm_delta exitflag])
|
||
% Visualisation
|
||
A = beta(1)*exp(-beta(2)*T);
|
||
subplot(2, 3, 3)
|
||
plot(T,A)
|
||
legend(txt_legend{1:i+2})
|
||
subplot(2, 3, 6)
|
||
plot(beta(1),beta(2),'ok')
|
||
text(beta(1),beta(2),[' \beta^{(' num2str(i) ')}'])
|
||
|
||
pause(0.5)
|
||
end
|
||
|
||
% Affichage des iteres de beta et sauvegarde des courbes
|
||
disp('--------------------------------------------------------------------------------------------')
|
||
print('C14_figures','-dpng')
|
||
diary
|