% This is the code for Problem 1 of HW 11.

% -------------------------------------------------------------------------
% The first part of the code is the same as code hw10_p1.m, except that the notations
%  have been changed to correspond to the notations of Lecture 11.
%  This part is used to produce a jagged function  uhat(t). 
% 

rand('state',10);   % This forces the random number generator to always use the same seed
                    %  (and hence produce the same output).
                   
tmax=20;              % length of time window
N=100;                % # of discretization points per this window, including its end points
h=tmax/(N-1);         % discretization step size
t=linspace(0,tmax,N);       % Vector of discretization points,
                            %  which has the same length as the function  Z  (see below).
                            
Z=log(t+5)+0.4*t;      % function Z without noise 
uexact=1./(t+5)+0.4;   % The exact  dZ/dt = u_exact -- this will be used as a benchmark for
                       %  the smoothened solution u(t) obtained in the 2nd part of the code.
                       
ep=0.005;                       % amplitude of noise
Zep = Z + ep*(1-2*rand(1,N));   % Z perturbed by a small noise of maximum magnitude  ep;
                                %  the jagged function  uhat  is derived from this noisy Z.
for k=1:N-1
	uhat(k)=(Zep(k+1)-Zep(k))/h;
end
uhat(N)=uhat(N-1);

% In this figure we will plot  uhat, uexact, and the smoothened solutions  u.
%  (The latter will be obtained in the second part of this code below.)
figure(1);
plot(t,uhat)
axis([0 tmax 0.4 0.62])
xlabel('t','fontsize',13)
hold on
plot(t,uexact,'g--','linewidth',3)

% -------------------------------------------------------------------------------
% This second part takes the jagged uhat(t) and smoothens it.
%

alpha=0.001;         % parameter  alpha  defined in the notes

der2Bmin=2;          % This allows us to choose which derivative is to be minimized:
                     % 2 => minimize the 2nd derivative, as in Lecture 11;
                     % 1 => minimize the 1st derivative (to compare with the above option).
                     
if  der2Bmin == 2
    
    % Set up matrix H defined in Eq. (16) of Lecture 11.
    %  If you want to see what this H looks like, change the N above to be some
    %  small number (within ten) and then display H by typing H at the Matlab's prompt.
    % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    H=diag(2*ones(1,N-2)) - diag(ones(1,N-3),1) - diag(ones(1,N-3),-1);
    H(1,:)=zeros(1,N-2);
    H(end,:)=zeros(1,N-2);
    H=H/h^2;
    M=transpose(H)*H;
    
    % Obtain the smoothened curve  u = (I+alpha*M)^(-1) * uhat:
    u = inv(eye(N-2)+alpha*M) * uhat(2:end-1)';    % the end points are treated next:
    u1=2*u(1)-u(2);
    uN=2*u(end)-u(end-1);
    u = [u1 u' uN];                   % augment the smoothemed  u  with the end points
    
    % Add the smoothened u curve to the curves of  uhat  and  uexact in the same figure.
    figure(1);
    plot(t,u,'r','linewidth',2)
    
    % Here we plot the differences between  uhat and uexact  and  u_smoothened and uexact.
    figure(2);
    xlabel('t','fontsize',13)
    plot(t,uhat-uexact)
    hold on
    plot(t,u-uexact,'r','linewidth',2)

    
elseif der2Bmin == 1
    % Set up matrix H corresponding to this case (not shown in the notes).
    % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    H=-diag(ones(1,N-1)) + diag(ones(1,N-2),1);
    H(end,end)=1;                 % this and the next lines account for condition u'_n=u'_(N-1).
    H(end,end-1)=-1;
    H=H/h;
    M=transpose(H)*H;
 
    % Obtain the smoothened curve  u = (I+alpha*M)^(-1) * uhat:
    u = inv(eye(N-1)+alpha*M) * uhat(1:end-1)';    % the right end point is treated next:
    uN=2*u(end)-u(end-1);
    u = [u' uN];                      % augment the smoothemed  u  with the right end point
    
    % Add the smoothened u curve to the curves of  uhat  and  uexact in the same figure.    
    figure(1);
    plot(t,u,'k','linewidth',2)

    % Here we plot the differences between  uhat and uexact  and  u_smoothened and uexact.
    figure(2);
    xlabel('t','fontsize',13)
    plot(t,uhat-uexact)
    hold on
    plot(t,u-uexact,'k','linewidth',2)
    
end

figure(1);
legend('uhat','uexact','usmoothned')
title([' \alpha = ' num2str(alpha)],'fontsize',13)
hold off
figure(2);
legend('uhat-uexact','usmoothned-uexact')
title([' \alpha = ' num2str(alpha)],'fontsize',13)
hold off


