Examples4

Tiedoston sisältö

% Jyvaskyl Summer School 2014
% Aku Seppanen
% University of Eastern Finland, Kuopio
% Department of Applied Physics

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 4.1: numerical differentiation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M = 20; % try with different discretizations: M = 5, 10, 20, 50, 100
x = linspace(0,2*pi,M+1)';
delta_x = x(2)-x(1);
y_true = sin(x);
q_true = cos(x);

var_n = .05; % variance of noise n
n = sqrt(var_n)*randn(length(x),1);
y = y_true + n;
figure(1),clf,subplot(2,1,1),plot(x,y_true,'g-',x,y,'+');
title('y=y(x)')

% Observation matrix
K = [ones(M+1,1),toeplitz([0; delta_x*(ones(M,1))],zeros(1,M))];

%%%  LS solution %%%
f_LS = inv(K'*K)*K'*y;
q_LS = f_LS(2:end);

x2 = x(1:end-1) + delta_x/2;
figure(1),subplot(2,1,2),plot(x,q_true,'g-',x2,q_LS,'r+');
title('q=q(x)')
legend('true q(x)','q_{LS}')

y_LS = K*f_LS;
figure(1),subplot(2,1,1),hold on, plot(x,y_LS,'ro');
legend('true y(x)','observed y','K*f_{LS}')


%%% Generalized Tikhonov regularization %%%
alpha = .5; % Try out different regularization parameter values
L_alpha = toeplitz(zeros(M-1,1),[0 -1 1 zeros(1,M-2)]); % 1st order difference matrix
%L_alpha = toeplitz(zeros(M-2,1),[0 1 -2 1 zeros(1,M-3)]); % 2nd order difference matrix
f_alpha = inv(K'*K + alpha^2*L_alpha'*L_alpha)*K'*y;
q_alpha = f_alpha(2:end);

figure(2),clf,subplot(2,1,1),plot(x,y_true,'g-',x,y,'+');
title('y=y(x)')
figure(2),subplot(2,1,2),plot(x,q_true,'g-',x2,q_alpha,'r+');
title(['q=q(x), regularized solution, \alpha = ',num2str(alpha)])
legend('true q(x)','q_{alpha}')

y_alpha = K*f_alpha;
figure(2),subplot(2,1,1),hold on, plot(x,y_alpha,'ro');
legend('true y(x)','observed y','K*f_{alpha}')