% 1D interpolation implemented with least-squares, varying epsilon

clear all;

Nx=101;
dx=1.0;
xmin=0.0;
xmax = xmin+dx*(Nx-1);
x = xmin+dx*[0:Nx-1]';

Nd=40;
rowindex=sort(unidrnd(Nx,Nd,1));
xd=xmin+dx*rowindex;
f=0.02;
y = sin(2*pi*f*xd);

figure(1);
clf;
plot(xd,y,'r+');
hold on;
axis([xmin xmax -2 2]);
xlabel('x');
ylabel('y');

M=Nx;
m=zeros(M,1);

N=Nd;
G=zeros(N,M);
d=zeros(N,1);
for p = [1:Nd]
    G(p,rowindex(p)) = 1.0;
    d(p)=y(p);
end

ipwr=0;
for pwr = [-6, -5, -4, -3, -2, -1, 0]
ipwr=ipwr+1;

epsilon = 10^pwr;
a=epsilon*(1.0);
b=epsilon*(-2.0);
c=epsilon*(1.0);
f=epsilon*(-1.0);

Nc=Nx;
D=zeros(Nc,M);
for p = [1:Nc-2]
    D(p,p) = a;
    D(p,p+1) = b;
    D(p,p+2) = c;
end

D(Nc-1,1)=a;
D(Nc-1,2)=f;

D(Nc,M-1)=a;
D(Nc,M)=f;

m= inv(G'*G + D'*D)*G'*d;
plot(x,m,'b');
hold on;

e = y-m(rowindex(:));
E(ipwr)=log10(e'*e);
end

figure(2);
clf;
plot([-6, -5, -4, -3, -2, -1, 0],E,'r');
hold on;
xlabel('log10 epsilon');
ylabel('log10 E');


