function y=tps(r)
% This is the thin-plate spline
if r < 0.000000000001
y=0;
else
y=r^2*log(r);
end
function y=fun(point)
% my target function
x=point(1);
z=point(2);
y=7-4*x^2+z^3;
function y=interpmat(points)
% this file computes the interpolation matrix
[n,m]=size(points);
for i=1:n
ipoint=points(i,:);
for j=1:n
jpoint=points(j,:);
mat(i,j)=tps(norm(ipoint-jpoint));
end
mat(i,n+1)=1;
mat(n+1,i)=1;
mat(i,n+2)=ipoint(1);
mat(n+2,i)=ipoint(1);
mat(i,n+3)=ipoint(2);
mat(n+3,i)=ipoint(2);
end
y=mat;
function y=righthandside(points)
% find the right hand side of the interpolation system
[n,m]=size(points);
for i=1:n
y(i)=fun(points(i,:));
end
y(n+1)=0;
y(n+2)=0;
y(n+3)=0;
function y=testinterp(n,m)
% this file test our interpolation routine
points=rand(n,2);
% compute interpolation matrix
mat=interpmat(points);
% create the right hand side
rhs=righthandside(points);
solution=inv(mat)*rhs';
% test solution
tpoints=rand(m,2)
maxerr=testsol(points,tpoints,solution);