Compute linear interpolant in pp form and find derivative

I would like to compute a linear interpolant in pp-form and then use the function fnder to compute its derivative. [The reason to do this is that I am solving a functional equation using collocation and I approximate the policy function using splines, but this is not so relevant here].

I used to do the following with interp1 and 'pp' (this is a minimum working example), but Matlab says that it will be removed in a later release. So my question is how to update this code.

Note: if I use cubic spline or other splines like pchip or makima I know how to do it, e.g.

pp = spline(x_grid,f_grid);

der_pp = fnder(pp,1);

y = ppval(der_pp,x)

but I want to do it for linear splines

Any help is greatly appreciated!

clear;clc;close all
n = 10;
x_grid = linspace(-5,5,n)';
x_grid_fine = linspace(-5,5,10*n)';
f_grid = x_grid.^2;
% Linear interpolant
pp = interp1(x_grid,f_grid,'linear','pp'); % MATLAB does not like it :(
% Derivative of linear interpolant
der_pp = fnder(pp,1);
fun_f_x = @(x) ppval(pp,x);
fun_df_dx = @(x) ppval(der_pp,x);
figure
plot(x_grid,fun_f_x(x_grid),'-o')
hold on
plot(x_grid_fine,fun_f_x(x_grid_fine))
legend('Raw data','Fitted')

figure
plot(x_grid,2*x_grid,'-o')
hold on
plot(x_grid_fine,fun_df_dx(x_grid_fine))
legend('Raw data','Fitted')

Answer 1

[John D'Errico]

n = 10;
x_grid = linspace(-5,5,n)';
f_grid = x_grid.^2; % A simple quadratic polynomial
fn = spapi(2,x_grid,f_grid)

fn = struct with fields:
      form: 'B-'
     knots: [-5 -5 -3.8889 -2.7778 -1.6667 -0.5556 0.5556 1.6667 2.7778 3.8889 5 5]
     coefs: [25 15.1235 7.7160 2.7778 0.3086 0.3086 2.7778 7.7160 15.1235 25]
    number: 10
     order: 2
       dim: 1

Convert to a pp-form. I might be able to do this directly in the call to spapi, but why bother?

pp = fn2fm(fn,'pp')

pp = struct with fields:
      form: 'pp'
    breaks: [-5 -3.8889 -2.7778 -1.6667 -0.5556 0.5556 1.6667 2.7778 3.8889 5]
     coefs: [9×2 double]
    pieces: 9
     order: 2
       dim: 1

fnplt(pp)

So clearly a piecewise linear spline. Now differentiate using fnder.

dfdx = fnder(pp)

dfdx = struct with fields:
      form: 'pp'
    breaks: [-5 -3.8889 -2.7778 -1.6667 -0.5556 0.5556 1.6667 2.7778 3.8889 5]
     coefs: [9×1 double]
    pieces: 9
     order: 1
       dim: 1

fnplt(dfdx)

So easy enough, especially if the interp1 usage is going to turn into a pumpkin one day.

[Alessandro Maria Marco]

Thank you so much for your answer!

I wanted to add a comment that might be useful for a general audience. Since the routine spapi returns the interpolant in B-form, one might wonder why I want instead the pp-form (which, after all, requires calling another routine for conversion, fn2fm). The reason is that the pp-form is much faster to evaluate.

This code

V_interp = spapi(2,k_grid,V0);
V_interp = fn2fm(V_interp,'pp');
% omitted
% - Evaluate linear spline in pp-form
V_interp_val = ppval(V_interp,kp_val);

is about 4 times faster on my laptop than this code that does the same:

V_interp = spapi(2,k_grid,V0);
% omitted
% - Evaluate linear spline in B-form
V_interp_val = fnval(V_interp,kp_val);