%% Bivariate Spline Interpolant and Osculatory Interpolation to Gridded Data
% The function spapi([0 0 0 0 1 2 2 2 2],[0 1 1 1 2],[2 0 1 2 -1]) produces the unique cubic spline f on the interval [0...2] with exactly one interior knot, at 1, that satisfies the five conditions
% f\left(0+\right)=2,f\left(1\right)=0,\textrm{Df}\left(1\right)=1,D^2 f\left(1\right)=2,f\left(2-\right)=-1\ldotp
% These include 3-fold matching at 1, i.e., matching there to prescribed values of the function and its first two derivatives.
% Here is an example of osculatory interpolation, to values y and slopes s at the sites x by a quintic spline:
% sp = spapi(augknt(x,6,2),[x,x,min(x),max(x)],[y,s,ddy0,ddy1]);
% with ddy0 and ddy1 values for the second derivative at the endpoints.
% As a related example, if you want to interpolate the sin(x) function at the distinct data sites by a cubic spline, and to match its slope at a subsequence x(s), then call the spapi function with these arguments:
% sp = spapi(4,[x x(s)], [sin(x) cos(x(s))]).
% The aptknt function will provide a suitable knot sequence. If you want to interpolate the same data by quintic splines, then simply change the value 4 to 6.
% As a bivariate example, here is a bivariate interpolant.

x = -2:.5:2; y = -1:.25:1; [xx, yy] = ndgrid(x,y);
z = exp(-(xx.^2+yy.^2));
sp = spapi({3,4},{x,y},z); fnplt(sp)

% As an illustration of osculatory interpolation to gridded data, here is complete bicubic interpolation, with the data explicitly derived from the bicubic polynomial g\left(u,v\right)=u^3 v^3. This is helpful to see exactly where the slopes, and slopes of slopes (the cross derivatives), must be placed in the data values supplied. Since g is a bicubic polynomial, its interpolant, f, must be g itself. Test this:

sites = {[0,1],[0,2]}; coefs = zeros(4,4); coefs(1,1) = 1;
g = ppmak(sites,coefs);
Dxg = fnval(fnder(g,[1,0]),sites);
Dyg = fnval(fnder(g,[0,1]),sites);
Dxyg = fnval(fnder(g,[1,1]),sites);
f = spapi({4,4}, {sites{1}([1,2,1,2]),sites{2}([1,2,1,2])}, ...
          [fnval(g,sites), Dyg ; ...
           Dxg.'         , Dxyg]);
if any( squeeze( fnbrk(fn2fm(f,'pp'), 'c') ) - coefs )
'something went wrong', end

% Copyright 2019 The MathWorks, Inc.