# Problem 4.6

meshfree 1D interpolation with compactly supported RBFs here: comparison of full vs. sparse version

f = @(x) 4*x.*(1-x);
% Wendlands \phi_3,2, full and sparse version
phi = @(r) max(1-r,0).^6.*(35*r.^2+18*r+3)/3;
phis = @(S,r) max(spones(S)-r,0).^6.*(35*r.^2+18*r+3*spones(S))/3;

N = 100;                       % vary this
X = linspace(0,1,N)'; y = f(X);
M = 2*N;
X_e = linspace(1/M,1-1/M,M)';   % has to be disjoint from X (sparsity structure)
ep = N/3;                       % vary this

% full version
tic
D = distm(X,X);
D = phi(ep*D);
c = D\y;
D_e = distm(X_e,X);
D_e = phi(ep*D_e);
y_e = D_e*c;
err = norm(y_e-f(X_e),inf)
plot(X_e,y_e,'o-'); hold on
toc

% sparse version
tic
D = sdistm(X,X,1/ep);
D = phis(D+speye(size(D)),ep*D);
c = D\y;
D_e = sdistm(X_e,X,1/ep);
D_e = phis(D_e,ep*D_e);
y_e = D_e*c;
err = norm(y_e-f(X_e),inf)
nnz(D)/size(D,1)                 % number of nonzeros per row in D
plot(X_e,y_e,'ro-'); hold off
toc

err =
0.0049
Elapsed time is 0.012996 seconds.
err =
0.0049
ans =
4.9400
Elapsed time is 0.023496 seconds. 