
% This code is a modified version of the code  hw10_p3.m  for 
%  Problems 3 of HW 10.
%  The modification concerns the part where one adds a smoothing matrix
%  alpha*B  to regularize the solution.

% ----------------------------------------------------------
% The first part of the code is essentially the same as in  hw10_p3.m.
% 

% (The length of the integration interval is scaled to equal 1 (one).)
n=20;            % # of discretization subintervals in the entire interval
h=1/n;           % discretization step size
xx=linspace(h/2,1-h/2,n);       % vector of discretization points,
                                %  located at the midpoints of the subintervals
                          
w=ones(1,n);  %(sin(pi*xx)).^2;     % gold density that engenders the exact solution mu(x) below

A=h./( (xx'*ones(1,n) - ones(n,1)*xx).^2 + 1 ).^(3/2);   
                            % matrix A of the discretized integral equation
mu=A*w';                    % the solution computed using the exact gold density

rand('state',0)
ep=0.0*max(abs(mu));       % amplitude of noise, set relative to the maximum 
muep=mu+ep*(1-2*rand(n,1)); % the solution perturbed by a small noise of maximum magnitude ep

                            
% ---------------------------------------------------------------
% This 2nd part is new. Regularization of the solution is performed here.
%

alpha=0.0;                % the regularization parameter  alpha, as in the notes

whichB=1;                   % This allows you to select a smoothing matrix.
                            % 1 => B = I (the identity matrix);
                            % 2 => B = Hmod defined in the notes
                            
if whichB == 1
    B=eye(n);
elseif whichB == 2
    B=diag(2*ones(1,n)) - diag(ones(1,n-1),1) - diag(ones(1,n-1),-1);
    B(1,1)=3;
    B(end,end)=3;
    B=B/h^2;
end
                 
% Compute the reconstructed solution:
wep=inv(A+alpha*B)*muep;           % gold density reconstructed from the solution  muep


figure(1);
plot(xx,wep,xx,w,'r--','linewidth',2)
legend('numerical','exact')
title('Numerical and exact gold densities w(x)','fontsize',16)

figure(2);
plot(xx,muep,xx,A*w','r--','linewidth',2)
legend('perturbed','exact')
title('Perturbed and unperturbed solutions mu(x)','fontsize',16)
