Skip to content
Snippets Groups Projects
Commit e01f1906 authored by Lisa Schätzle's avatar Lisa Schätzle
Browse files

Update example_5_2.m

parent 002e5d92
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 75 additions and 52 deletions
src/example_5_2.m
+ 75
− 52
View file @ e01f1906
...@@ -18,7 +18,7 @@ q = [-.5, 1]; ...@@ -18,7 +18,7 @@ q = [-.5, 1];
noiselevel = 0.05; 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); nr_rep = length(alpha_vec);
...@@ -26,7 +26,7 @@ P_Omega = cell(1,nr_rep); ...@@ -26,7 +26,7 @@ P_Omega = cell(1,nr_rep);
P_OmegaC = cell(1,nr_rep); P_OmegaC = cell(1,nr_rep);
for iteri=1:nr_rep for iteri=1:nr_rep
alpha = alpha_vec(iteri); alpha = alpha_vec(iteri);
% cross-shaped Omega, 20% midding data: % cross-shaped Omega, 20% midding data:
...@@ -38,6 +38,9 @@ for iteri=1:nr_rep ...@@ -38,6 +38,9 @@ for iteri=1:nr_rep
% To enforce symmetry according to reciprocity principle: % To enforce symmetry according to reciprocity principle:
P_Omega{iteri} = applyRP(Help); P_Omega{iteri} = applyRP(Help);
P_OmegaC{iteri} = ones(sampling.nd,sampling.nxhat) - P_Omega{iteri}; P_OmegaC{iteri} = ones(sampling.nd,sampling.nxhat) - P_Omega{iteri};
ratio(iteri) = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
clear Help clear Help
end end
...@@ -64,28 +67,21 @@ end ...@@ -64,28 +67,21 @@ end
% Numerical reconstruction with least squares approach: % Numerical reconstruction with least squares approach:
tol = 1e-5; tol = [1e-5 1e-5 1e-5 1e-4 1e-4 5e-3 5e-3 5e-3 1e-3 1e-3];
kmax = 100; kmax = 300;
Aall = cell(1,nr_rep); Aall = cell(1,nr_rep);
B = cell(1,nr_rep); B = cell(1,nr_rep);
for iteri = 1:nr_rep 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}; Aall{iteri} = A{2};
B{iteri} = A{1}; B{iteri} = A{1};
clear A clear A
relerrOmega(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall); relerrOmega1_cg(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall);
relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); relerr1_cg(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
end
disp('Table 5.2 (left): Relative errors of far field operator completion 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))
end end
% Numerical reconstruction with l1xl1-minimization: % Numerical reconstruction with l1xl1-minimization:
...@@ -96,22 +92,37 @@ Aall = cell(1,nr_rep); ...@@ -96,22 +92,37 @@ Aall = cell(1,nr_rep);
B = cell(1,nr_rep); B = cell(1,nr_rep);
for iteri = 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}; Aall{iteri} = A{1};
clear A clear A
relerrOmega(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall); relerrOmega1_fista(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall);
relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); relerr1_fista(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
end end
disp('Table 5.2 (left): Relative errors of far field operator completion with l1xl1-minimization for Omega_1:') %% Plot of relative errors for Omega_1:
fprintf(' |Omega|/(4pi^2) epsilon epsilon_Omega\n')
for iteri = 1:nr_rep figure()
ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
fprintf(' %.2f %.3f %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3)) semilogy(ratio, relerrOmega1_cg, 'r--o', ratio, relerr1_cg, 'b--*', ratio, relerrOmega1_fista, 'r-o', ratio, relerr1_fista, 'b-*', 'LineWidth', 1.2)
end
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')
% Plot of geometry: 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_completion.eps -depsc
%% Plot of geometry:
figure() figure()
...@@ -146,12 +157,15 @@ title('Geometry and a priori information', 'Interpreter', 'Latex') ...@@ -146,12 +157,15 @@ title('Geometry and a priori information', 'Interpreter', 'Latex')
grid on grid on
axis equal axis equal
text(26.2,-1.8,'$D_1$','Color','blue','FontSize',18, 'Interpreter','latex')
text(6.25,-10.8,'$D_2$','Color','blue','FontSize',18, 'Interpreter','latex')
ax = gca; ax = gca;
ax.FontSize = 17; ax.FontSize = 18;
print ../figures/geometry_vary_Omega.eps -depsc print ../figures/geometry_vary_Omega.eps -depsc
% Plot of Omega_1: %% Plot of Omega_1:
figure() figure()
...@@ -204,8 +218,6 @@ print ../figures/support_Omega1.eps -depsc ...@@ -204,8 +218,6 @@ print ../figures/support_Omega1.eps -depsc
%% Example 5.2 (Completion only) for Omega_2: %% Example 5.2 (Completion only) for Omega_2:
clear all
% Parameters: % Parameters:
k = .5; % wave number k = .5; % wave number
...@@ -220,7 +232,7 @@ q = [-.5, 1]; ...@@ -220,7 +232,7 @@ q = [-.5, 1];
noiselevel = 0.05; 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); nr_rep = length(alpha_vec);
...@@ -240,6 +252,9 @@ for iteri=1:nr_rep ...@@ -240,6 +252,9 @@ for iteri=1:nr_rep
% To enforce symmetry according to reciprocity principle: % To enforce symmetry according to reciprocity principle:
P_Omega{iteri} = applyRP(Help); P_Omega{iteri} = applyRP(Help);
P_OmegaC{iteri} = ones(sampling.nd,sampling.nxhat) - P_Omega{iteri}; P_OmegaC{iteri} = ones(sampling.nd,sampling.nxhat) - P_Omega{iteri};
ratio(iteri) = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
clear Help clear Help
end end
...@@ -266,28 +281,21 @@ end ...@@ -266,28 +281,21 @@ end
% Numerical reconstruction with least squares approach: % Numerical reconstruction with least squares approach:
tol = 1e-5; tol = [1e-5 1e-5 1e-5 1e-4 5e-3 5e-3 5e-3 5e-3 1e-3 1e-3];
kmax = 100; kmax = 300;
Aall = cell(1,nr_rep); Aall = cell(1,nr_rep);
B = cell(1,nr_rep); B = cell(1,nr_rep);
for iteri = 1:nr_rep 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}; Aall{iteri} = A{2};
B{iteri} = A{1}; B{iteri} = A{1};
clear A clear A
relerrOmega(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall); relerrOmega2_cg(iteri) = norm(B{iteri}+P_Omega{iteri}.*Fall)/norm(P_Omega{iteri}.*Fall);
relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); relerr2_cg(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
end
disp('Table 5.2 (right): Relative errors of far field operator completion 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))
end end
% Numerical reconstruction with l1xl1-minimization: % Numerical reconstruction with l1xl1-minimization:
...@@ -298,22 +306,37 @@ Aall = cell(1,nr_rep); ...@@ -298,22 +306,37 @@ Aall = cell(1,nr_rep);
B = cell(1,nr_rep); B = cell(1,nr_rep);
for iteri = 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}; Aall{iteri} = A{1};
clear A clear A
relerrOmega(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall); relerrOmega2_fista(iteri) = norm(P_Omega{iteri}.*(Aall{iteri}-Fall))/norm(P_Omega{iteri}.*Fall);
relerr(iteri) = norm(Aall{iteri}-Fall)/norm(Fall); relerr2_fista(iteri) = norm(Aall{iteri}-Fall)/norm(Fall);
end end
disp('Table 5.2 (right): Relative errors of far field operator completion with l1xl1-minimization for Omega_2:') %% Plot of relative errors for Omega_2:
fprintf(' |Omega|/(4pi^2) epsilon epsilon_Omega\n')
for iteri = 1:nr_rep figure()
ratio = sum(sum(P_Omega{iteri}))/(sampling.nxhat*sampling.nd);
fprintf(' %.2f %.3f %.3f\n', ratio, round(relerr(iteri),3), round(relerrOmega(iteri),3)) semilogy(ratio, relerrOmega2_cg, 'r--o', ratio, relerr2_cg, 'b--*', ratio, relerrOmega2_fista, 'r-o', ratio, relerr2_fista, 'b-*', 'LineWidth', 1.2)
end
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)])
ylim([10^(-2) 1])
grid on
ax = gca;
ax.FontSize = 16;
print ../figures/errors_vary_Omega2_completion.eps -depsc
% Plot of Omega_2: %% Plot of Omega_2:
figure() figure()
...@@ -385,4 +408,4 @@ title('Non-observable set $\Omega_2$', 'Interpreter', 'Latex') ...@@ -385,4 +408,4 @@ title('Non-observable set $\Omega_2$', 'Interpreter', 'Latex')
set(gca,'Fontsize',18) set(gca,'Fontsize',18)
set(gca,'YDir','normal') set(gca,'YDir','normal')
print ../figures/support_Omega2.eps -depsc print ../figures/support_Omega2.eps -depsc
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Help | Imprint | Privacy policy | Accessibility | Contact