Commit e7c52e2a authored by Gidon Lucian Bauer's avatar Gidon Lucian Bauer 🏳
Browse files

MSM Bonus 03 fertig bis auf schöne Ausgabe

parent 1587805e
......@@ -33,9 +33,9 @@ load error.mat; % error (random noise)
u_f_meas = u_f + error; % measurement of final temperature
u_0_est = expm(T*C)*u_f_meas;
sprintf('1d backward heat equation with n=%i and noise level sigma=%g\n==============================================================', n, sigma)
sprintf('1d backward heat equation with n=%i and noise level sigma=%g\n==============================================================', n, sigma);
err = norm(u_0 - u_0_est)
err = norm(u_0 - u_0_est);
figure(1);
subplot(1,2,1);
......@@ -58,7 +58,7 @@ if (beta <= 0 || beta >= s(1)^-2)
warning('Beta hat einen falschen Wert.'); % check if beta has a correct value (a)
end
k_star = floor(max_iter/2); % stopping index
k_star = -1; % stopping index
eps_alpha = 1e-01;
disp(d(k_star, A, u_f_meas));
......@@ -71,20 +71,36 @@ disp(d(k_star, A, u_f_meas));
% k_star = i;
% end
for i = 2:min(k_star, max_iter)
xi_norm = zeros(max_iter-1, 1);
disc = zeros(max_iter-1, 1);
C_xi_norm = zeros(max_iter-1, 1);
for i = 2:max_iter
iter(:,i) = iter(:,i-1) - beta * A' * (A * iter(:,i-1) - u_f_meas);
% c)?
xi_norm(i-1) = norm(iter(:,i), 2);
disc(i-1) = norm(A*iter(:,i) - u_f_meas, 2);
% d)
C_xi_norm(i-1) = norm(C*iter(:,i), 2);
% b)
if disc(i-1) <= eps_alpha && k_star == -1
k_star = i;
disp(k_star);
disp(norm(A*iter(:,i) - u_f_meas, 2));
end
end
% b)
figure(2);
plot(x, u_0, 'DisplayName', 'exact');
hold on;
plot(x, iter(:, k_star), 'DisplayName', 'estimate');
plot(x, iter(:, k_star), 'DisplayName', 'estimate for k*');
plot(x, iter(:, end), 'DisplayName', 'estimate max iter');
xlabel('x');
ylabel('temperature');
legend;
hold off;
% b)
figure(3);
plot(x, u_f, 'DisplayName', 'exact');
hold on;
......@@ -94,3 +110,23 @@ xlabel('x');
ylabel('temperature');
legend;
hold off;
figure(4);
% loglog(2:max_iter, disc, '*', 'DisplayName', 'Discrepancy');
loglog(disc(1:end), xi_norm(1:end), '*', 'DisplayName', 'Discrepancy');
hold on;
% loglog(2:max_iter, temp, 'DisplayName', '\xi');
xlabel('Discrepancy');
ylabel('\xi norm');
title('L curve');
legend;
hold off;
figure(5);
loglog(disc, C_xi_norm, '*', 'DisplayName', 'Discrepancy');
hold on;
xlabel('Discrepancy');
ylabel('C-\xi norm');
title('L curve');
legend;
hold off;
\ No newline at end of file
function res = d_der(alpha, s, y_tilde)
end
\ No newline at end of file
Markdown is supported
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