% 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}')