diff --git a/src/example_5_3.m b/src/example_5_3.m index 5d1a8cc871cb33c08ae6968b3a59844e13510240..08ef90f8b0f8af29ce2e973a81055b841d16a6d4 100644 --- a/src/example_5_3.m +++ b/src/example_5_3.m @@ -18,7 +18,7 @@ q = [-.5, 1]; noiselevel = 0.05; -alpha_vec = pi/16*[1,2,3,4,5,6]; % controlls area of non-observable sets +alpha_vec = pi/32*[1,2,3,4,5,6,7,8,9,10]; % controlls area of non-observable sets nr_rep = length(alpha_vec); @@ -38,6 +38,9 @@ for iteri=1:nr_rep % To enforce symmetry according to reciprocity principle: P_Omega{iteri} = applyRP(Help); P_OmegaC{iteri} = ones(sampling.nd,sampling.nxhat) - P_Omega{iteri}; + + ratio(iteri) = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd); + clear Help end @@ -64,8 +67,8 @@ end % Numerical reconstruction with least squares approach: -tol = 1e-5; -kmax = 100; +tol = [1e-5 1e-5 1e-5 1e-4 1e-4 5e-3 5e-3 5e-3 5e-3 1e-3]; +kmax = 300; Anut = cell(1,nr_rep); Akite = cell(1,nr_rep); @@ -75,7 +78,7 @@ B = cell(1,nr_rep); indOfInterest = sub2ind([nr_of_scatterers,nr_of_scatterers],1:nr_of_scatterers,1:nr_of_scatterers); for iteri = 1:nr_rep - [A, ~, ~] = CG_secondOrder(G{iteri}, sampling, k, P_Omega{iteri}, z, R, kmax, tol); + [A, ~, ~] = CG_secondOrder(G{iteri}, sampling, k, P_Omega{iteri}, z, R, kmax, tol(iteri)); Aall{iteri} = A{2} + A{3} + A{4} + A{5}; A = {A{1},A{indOfInterest+1}}; @@ -85,15 +88,8 @@ for iteri = 1:nr_rep clear A - relerrOmega(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall); - relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); -end - -disp('Table 5.3 (left): Relative errors of far field operator completion and splitting with least squares approach for Omega_1:') -fprintf(' |Omega|/(4pi^2) epsilon epsilon_Omega\n') -for iteri = 1:nr_rep - ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd); - fprintf(' %.2f %.3f %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3)) + relerrOmega1_cg(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall); + relerr1_cg(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); end % Numerical reconstruction with l1xl1-minimization: @@ -106,23 +102,38 @@ Aall = cell(1,nr_rep); B = cell(1,nr_rep); for iteri = 1:nr_rep - [A, nr_of_iterations, Resnorm] = FISTA_secondOrder(G{iteri}, sampling, k, P_OmegaC{iteri}, z, kmax); + [A, ~, ~] = FISTA_secondOrder(G{iteri}, sampling, k, P_OmegaC{iteri}, z, kmax); Aall{iteri} = A{1}; Anut{iteri} = A{2}; Akite{iteri} = A{3}; clear A - relerrOmega(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall); - relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); + relerrOmega1_fista(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall); + relerr1_fista(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); end -disp('Table 5.3 (left): Relative errors of far field operator completion and splitting with l1xl1-minimization for Omega_1:') -fprintf(' |Omega|/(4pi^2) epsilon epsilon_Omega\n') -for iteri = 1:nr_rep - ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd); - fprintf(' %.2f %.3f %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3)) -end +%% Plot of relative errors for Omega_1: + +figure() + +semilogy(ratio, relerrOmega1_cg, 'r--o', ratio, relerr1_cg, 'b--*', ratio, relerrOmega1_fista, 'r-o', ratio, relerr1_fista, 'b-*', 'LineWidth', 1.2) + +title('Relative errors for varying area of non-observable set $\Omega_1$', 'Interpreter', 'latex') +xlabel('$\frac{|\Omega|}{4\pi^2}$', 'Interpreter', 'latex') +ylabel('$\epsilon_{\mathrm{rel}}$', 'Interpreter', 'latex') + +legend({'$\epsilon_{\mathrm{rel}}^{\Omega}$ using cg','$\epsilon_{\mathrm{rel}}$ using cg', '$\epsilon_{\mathrm{rel}}^{\Omega}$ using FISTA', '$\epsilon_{\mathrm{rel}}$ using FISTA'}, 'Interpreter','latex', 'Location','southeast') + +xlim([ratio(1) ratio(end)]) +ylim([10^(-2) 1]) + +grid on + +ax = gca; +ax.FontSize = 16; + +print ../figures/errors_vary_Omega1_completionSplitting.eps -depsc %% Example 5.3 (Completion and Splitting) for Omega_2: @@ -142,7 +153,7 @@ q = [-.5, 1]; noiselevel = 0.05; -alpha_vec = pi/16*[1,2,3,4,5,6]; % controlls area of non-observable sets +alpha_vec = pi/32*[1,2,3,4,5,6,7,8,9,10]; % controlls area of non-observable sets nr_rep = length(alpha_vec); @@ -162,6 +173,9 @@ for iteri=1:nr_rep % To enforce symmetry according to reciprocity principle: P_Omega{iteri} = applyRP(Help); P_OmegaC{iteri} = ones(sampling.nd,sampling.nxhat) - P_Omega{iteri}; + + ratio(iteri) = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd); + clear Help end @@ -188,8 +202,8 @@ end % Numerical reconstruction with least squares approach: -tol = 1e-5; -kmax = 100; +tol = [1e-5 1e-5 1e-5 1e-4 1e-4 5e-3 5e-3 1e-3 1e-3 1e-3]; +kmax = 300; Anut = cell(1,nr_rep); Akite = cell(1,nr_rep); @@ -199,7 +213,7 @@ B = cell(1,nr_rep); indOfInterest = sub2ind([nr_of_scatterers,nr_of_scatterers],1:nr_of_scatterers,1:nr_of_scatterers); for iteri = 1:nr_rep - [A, ~, ~] = CG_secondOrder(G{iteri}, sampling, k, P_Omega{iteri}, z, R, kmax, tol); + [A, ~, ~] = CG_secondOrder(G{iteri}, sampling, k, P_Omega{iteri}, z, R, kmax, tol(iteri)); Aall{iteri} = A{2} + A{3} + A{4} + A{5}; A = {A{1},A{indOfInterest+1}}; @@ -209,15 +223,8 @@ for iteri = 1:nr_rep clear A - relerrOmega(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall); - relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); -end - -disp('Table 5.3 (right): Relative errors of far field operator completion and splitting with least squares approach for Omega_2:') -fprintf(' |Omega|/(4pi^2) epsilon epsilon_Omega\n') -for iteri = 1:nr_rep - ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd); - fprintf(' %.2f %.3f %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3)) + relerrOmega2_cg(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall); + relerr2_cg(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); end % Numerical reconstruction with l1xl1-minimization: @@ -230,20 +237,34 @@ Aall = cell(1,nr_rep); B = cell(1,nr_rep); for iteri = 1:nr_rep - [A, nr_of_iterations, Resnorm] = FISTA_secondOrder(G{iteri}, sampling, k, P_OmegaC{iteri}, z, kmax); + [A, ~, ~] = FISTA_secondOrder(G{iteri}, sampling, k, P_OmegaC{iteri}, z, kmax); Aall{iteri} = A{1}; Anut{iteri} = A{2}; Akite{iteri} = A{3}; clear A - relerrOmega(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall); - relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); + relerrOmega2_fista(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall); + relerr2_fista(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); end -disp('Table 5.3 (right): Relative errors of far field operator completion and splitting with l1xl1-minimization for Omega_2:') -fprintf(' |Omega|/(4pi^2) epsilon epsilon_Omega\n') -for iteri = 1:nr_rep - ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd); - fprintf(' %.2f %.3f %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3)) -end \ No newline at end of file +%% Plot of relative errors for Omega_2: + +figure() + +semilogy(ratio, relerrOmega2_cg, 'r--o', ratio, relerr2_cg, 'b--*', ratio, relerrOmega2_fista, 'r-o', ratio, relerr2_fista, 'b-*', 'LineWidth', 1.2) + +title('Relative errors for varying area of non-observable set $\Omega_2$', 'Interpreter', 'latex') +xlabel('$\frac{|\Omega|}{4\pi^2}$', 'Interpreter', 'latex') +ylabel('$\epsilon_{\mathrm{rel}}$', 'Interpreter', 'latex') + +legend({'$\epsilon_{\mathrm{rel}}^{\Omega}$ using cg','$\epsilon_{\mathrm{rel}}$ using cg', '$\epsilon_{\mathrm{rel}}^{\Omega}$ using FISTA', '$\epsilon_{\mathrm{rel}}$ using FISTA'}, 'Interpreter','latex', 'Location','southeast') + +xlim([ratio(1) ratio(end)]) + +grid on + +ax = gca; +ax.FontSize = 16; + +print ../figures/errors_vary_Omega2_completionSplitting.eps -depsc \ No newline at end of file