Cubic interpolation coefficients and basis matrix

Suppose I have n data points (x(i),y(i)) with i=1,..,n. I want to compute a cubic interpolant that fits exactly these points (interpolation, not least squares fit). I can do the following in Matlab:

x_min = -1;
x_max = 1;
nx    = 10;
x = linspace(x_min,x_max,nx)';
y = exp(-x);
pp = csapi(x,y); %Can also use spline
disp(pp)

      form: 'pp'
    breaks: [-1 -0.7778 -0.5556 -0.3333 -0.1111 0.1111 0.3333 0.5556 0.7778 1]
     coefs: [9×4 double]
    pieces: 9
     order: 4
       dim: 1

This returns the coefficients in the matrix pp.coefs: each row l (for l=1,..,n-1) of this matrix gives the 4 coefficients of the cubic polynomial for the specific subinterval l. However, I would like the cubic in another form. A cubic spline can be written as

(1)

where phi_j(x) are the basis functions and c is a vector of n+2 coefficients. How can I get these n+2 coefficients (and, optionally the basis matrix Phi)? The Matlab functions csapi and spline give this (n-1)*4 matrix of coefficients which is not what I want.

Reference: The source for (1) is Fehr and Kindermann, "Computational Economics", Oxford University Press, page 93.

Any help is greatly appreciated, thanks!

[Alessandro Maria Marco]

Please let me know in case the question is not clear or is not well-posed and if further clarifications are needed

[Torsten]

I don't feel like googling the equational form of the basis functions phi_j. Can you share a link ?

[Alessandro Maria Marco]

Sure, here is a wiki link: https://en.wikipedia.org/wiki/B-spline

I also attach the relevant part: any spline function can be written as a linear combination of basis functions:

[Alessandro Maria Marco]

The basis functions are defined recursively. This stuff is known as B-splines, I think Matlab should have it

[Torsten]

Then you will have to sum the B_(j,3) functions that have non-zero support in [t_i,t_i+1] and equate this sum to the usual cubic spline representation coming from "csape" in [t_i,t_i+1]. Comparison of coefficients of the two cubic polynomials in all subintervals should give you a linear system of equations to determine the alpha_i coeffcients from the "csape" coefficients.

Answer 1

[Matt J]

Using this FEX downloadable,

https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form

x_min = -1;
x_max = 1;
nx    = 10;
x = linspace(x_min,x_max,nx)';
xx=linspace(x(1),x(end),1000)';


fun=@(in)csapi(x,in,xx);

Basis=func2mat(fun,x); %Columns of this matrix are basis functions

y=exp(-x);
yy=exp(-xx);
c=Basis\yy; %Coefficients

subplot(1,2,1)
h=plot(x,y,'o',xx,Basis*c,'.r'); axis square
title({'Interp by Basis'; 'Matrix Multiplication'})

subplot(1,2,2)
plot(Basis);  axis square
title 'Basis functions'

[Alessandro Maria Marco]

Thanks for your answer! I wanted to note that there is another way to do this using the routine "spcol" from the curve fitting toolbox. Might be useful for those users who have the curve fitting toolbox and prefer not to download extra functions from FEX.