207 lines
7.2 KiB
Mathematica
207 lines
7.2 KiB
Mathematica
|
%-------------------------------------------------------------------------%
|
|||
|
% 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
|