clear all;
format long
%%%%%%%%%%%%%%%%%%%%%%%%Exercise 5.2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%collocation from exercise 3.1
%part a

%y and d=d(y) are data from Table 3.1 in text
y=[.025 .075 .125 .175 .225 ...
    .2750 .3250 .375 .425 .4725 ...
    .525 .575 .625 .675 .7225 ...
    .775 .825 .875 .925 .975]';
d=[.2388 .2319 .2252 .2188 .2126 ...
    .2066 .2008 .1952 .1898 .1846 ...
    .1795 .1746 .1699 .1654 .161 ...
    .1567 .1526 .1486 .1447 .141]';

%m is number of data points
%n is # of collocation points and a\leq x\leq b
%xi is the midpoint
n=20; %as n increases, condition number blows up, and...
     %thus our data may be ill-conditioned.
     %However, small values of n yield "not so accurate" norms
m=20; a=0; b=1; x=[];
deltax=(b-a)/n;

%generate G matrix so we can solve for m(x)
for i=1:m
    for j=1:n
        x=[x,a+(j-.5)*(deltax)];
        G(i,j)=(x(j)*exp(-x(j)*y(i)))*deltax;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%5.2 (e) Second-order Tikhonov%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

L2 = get_l(n,2);
[U2,sm2,X2,V2]=cgsvd(G,L2);

%GCV Method%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
alpha_secondgcv=gcv(U2,sm2,d);
alpha_secondgcv % alpha for GCV

[mtik_secondgcv,rho_secondgcv,eta_secondgcv]=tikhonov(U2,sm2,X2,d,alpha_secondgcv);
%mtik_secondgcv

%Discrepancy Principle%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use the discrepancy principle to get a second solution.

delta=(10^(-4))*norm(d);
[mtik_seconddisc,alpha_seconddisc]=discrep(U2,sm2,X2,d,delta);
alpha_seconddisc %alpha for discrepancy

[mtik_seconddisc,rho_seconddisc,eta_seconddisc]=tikhonov(U2,sm2,X2,d,alpha_seconddisc);
%mtik_seconddisc

%L-Curve Criterion%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First, plot the l_curve and find its corner.

figure
[alpha_secondL,rho_secondL,eta_secondL]=l_curve(U2,sm2,d);
alpha_secondL %alpha for L curve

[mtik_secondL,rho_secondL,eta_secondL]=tikhonov(U2,sm2,X2,d,alpha_secondL);
%mtik_secondL

%Plot all 3 solutions together with data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure
plot(mtik_secondgcv,'r')
hold on
plot(mtik_secondL,'k')
legend('mtik_ secondgcv','mtik_ secondL')
hold off

figure
plot(mtik_seconddisc,'b')
legend('mtik_ seconddisc')