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

Modellgestützte Schätzmethoden Bonus 03

parent e607c54e
% Backward heat equation in 1D
% Bonus problem 3: iterative regularization + parameter choice
close all;
clear;
L = pi; % length of the rod
T = 5; % final time
sigma = 1e-2; % magnitude of error
kappa = 1e-3; % heat conduction coefficient
n = 100;
N = n+1; % number of intervals
h = L/N; % grid size
% create C
C = zeros(n,n);
% first and last row
C(1,1:2) = [2 -1];
C(n,n-1:n) = [-1 2];
% inner rows
for i = 2:n-1
C(i,i-1:i+1) = [-1 2 -1];
end
% and, finally, multiply with the right factor
C = C*(kappa/h^2);
%% backward heat conduction
x = linspace(h,L-h,n)'; % generate linearly spaced vector
u_0 = sin(x) + 0.5*sin(4*x); % initial temperature
A = expm(-T*C);
u_f = A*u_0; % final temperature
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)
err = norm(u_0 - u_0_est)
figure(1);
subplot(1,2,1);
plot(x, u_f, x, u_f_meas, x, A*u_0_est);
legend('u_f', 'u_f measured', 'u_f est. (data fit)');
subplot(1,2,2);
plot(x, u_0, x, u_0_est, 'r');
legend('u_0', 'estimation for u_0');
% Landweber iterations
max_iter = 250;
iter = zeros(n,max_iter); % store all iterates
% ...
s = svd(A); % singular values of A
beta = 1;
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
eps_alpha = 1e-01;
disp(d(k_star, A, u_f_meas));
% d(alpha)
% d'(alpha)
% -> Eulermethode
% disc_para = norm(A * iter(:,i) - u_f_meas, 2);
% if (abs(disc_para - eps_alpha) <= 1e-8)
% k_star = i;
% end
for i = 2:min(k_star, max_iter)
iter(:,i) = iter(:,i-1) - beta * A' * (A * iter(:,i-1) - u_f_meas);
end
figure(2);
plot(x, u_0, 'DisplayName', 'exact');
hold on;
plot(x, iter(:, k_star), 'DisplayName', 'estimate');
xlabel('x');
ylabel('temperature');
legend;
hold off;
figure(3);
plot(x, u_f, 'DisplayName', 'exact');
hold on;
plot(x, u_f_meas, 'DisplayName', 'measurment');
plot(x, A*iter(:, k_star), 'DisplayName', 'y bar');
xlabel('x');
ylabel('temperature');
legend;
hold off;
function res = d(alpha, A, y_tilde)
% to_min = @(xi) norm(A * xi - y_tilde, 2)^2 + alpha^2 * norm(xi, 2)^2;
% xi_alpha = fminsearch(to_min, zeros(length(A)));
xi_alpha = (A'*y_tilde)\(A'*A + alpha^2*eye(size(A)));
xi_alpha = xi_alpha';
res = norm(A * xi_alpha - y_tilde, 2);
end
\ No newline at end of file
function res = d_der(alpha, s, y_tilde)
end
\ No newline at end of file
error = [0.0047586058727834262; 0.014122326864445039; 0.00022608484309598147; ...
-0.00047869410220205893; 0.01701334654274959; -0.0050971171276742652; ...
-2.8549601443212068E-5; 0.0091986707980639543; 0.0014980873263276118; ...
0.014049334456769766; 0.010341215395697101; 0.0029157028877080573; ...
-0.0077769853831160338; 0.0056669609753930471; -0.013826211594803522; ...
0.0024447467558988776; 0.0080843880316769125; 0.0021304169841700046; ...
0.0087967716423348885; 0.020388762514140422; 0.00923932448688937; ...
0.00266917446595828; 0.0064166150642177931; 0.0042548535562529554; ...
-0.013147235034868837; -0.0041641121969943386; 0.012246878247853367; ...
-0.0004358420554633264; 0.0058242327744796917; -0.010065000746193356; ...
0.00064516742311104412; 0.0060029194918578355; -0.013615149548645669; ...
0.0034759263196006453; -0.0018184321845933357; -0.0093953476594149151; ...
-0.00037533188819171594; -0.0189630449362245; -0.021279767681826739; ...
-0.011769233307149575; -0.009905322204841761; -0.011730323272674051; ...
-0.017254277895286919; 0.0028822808966501058; -0.015941837202668067; ...
0.0011021884922338119; 0.0078706667635798019; -2.2267863138361347E-5; ...
0.00093108760058803529; -0.0037815705658975751; -0.014826761110590034; ...
-0.000438185853582948; 0.0096082521168211459; 0.017382449326133406; ...
-0.0043020624242610045; -0.0162732273646978; 0.0016634749246006552; ...
0.0037626591045071883; -0.0022695046470623255; -0.0114891228961879; ...
0.0202433259954387; -0.023595235109584913; -0.0050997204579993904; ...
-0.013216255908101016; -0.0063612824966004147; 0.0031785141905969674; ...
0.0013804797444952586; -0.00710735074811226; 0.0077700352671978791; ...
0.0062239392417201287; 0.0064738088451604678; -0.0042563168166035135; ...
0.010485807605364415; 0.0066070708636717492; 0.025087724731851058; ...
0.01063459639041023; 0.011569216533273936; 0.00052978827325612845; ...
-0.012883860029402915; -0.0037122124749763225; -0.0075779191649881861; ...
-0.0056396891702742685; 0.0055513856000184636; -0.0055677806395036571; ...
-0.0089511313548243567; -0.0040932772163407042; -0.0016088676673844289; ...
0.0040933443045545531; -0.009526359971199775; 0.0031731747232461188; ...
0.0007802008097807753; 0.013243854491602236; -0.0021317048642158594; ...
-0.001344786445554264; -0.01171355817920296; -0.013852626831952397; ...
0.0031050831860049873; -0.0024948906447393891; 0.0050374405573732819; ...
-0.0089266140347703787];
function y = g_alpha(z, alpha)
y = z^2 / (z^2 * alpha^2);
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