MATLAB Answers

Bhomik Luthra

Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) for given data and then finding area

Asked by Bhomik Luthra
on 17 Jun 2013

Hello all

    I am trying to do Piecewise Cubic Hermite Interpolation on the data given below and then I want to get the area covered by the polynomials with x axis. I think, I am misunderstanding the meaning of coefficients returned by pchip command, but not sure. Does anyonw know what could be the problem? 

x = [5.8808 6.5137 7.1828 7.8953];

y = [31.2472 33.9977 36.7661 39.3567];

pp = pchip(x,y)

If I see see pp,it gives pp as

   form: 'pp'
    breaks: [5.8808 6.5137 7.1828 7.8953]
     coefs: [3x4 double]
    pieces: 3
     order: 4
       dim: 1

and pp.coefs are

-0.0112 -0.1529 4.4472 31.2472

-0.3613 0.0884 4.2401 33.9977

-0.0422 -0.3028 3.8731 36.7661

I think these are the polynomials representing the three intervals




But when I find the y values corresponding to x values using these polynomials, it gives wrong values.

It gives negative y values for second polynomial . Even third polynomial do not seem to satisfy the points.

I used these commands for obtaining the values

For example:- (for second polynomial)

xs = linspace(6.5137, 7.1828, 200);

y = polyval(pp.coefs(2,:),xs);


I want to find the area under the curve covered by this plot, that's why I am trying to find the polynomial. Is there any other way to do it or if anyone could find the problem in the commands that I am using, please let me know.


Bhomik Luthra


1 Answer

Answer by Andrei Bobrov
on 17 Jun 2013
Edited by Andrei Bobrov
on 17 Jun 2013
 Accepted answer
x = [5.8808 6.5137 7.1828 7.8953];
y = [31.2472 33.9977 36.7661 39.3567];
pch = pchip(x,y);
out = fnval(fnint(pch),x([2,3]))*[-1;1];    % if you have 'Curve Fitting Toolbox'


out = integral(@(x)ppval(pch,x),x(2),x(3)); % 


Bhomik Luthra
on 17 Jun 2013

Thanks for your valuable answer, but I still have the same doubt. What exactly are pch.coefs? Could you please tell me what do they represent?

Thanks Bhomik

Andrei Bobrov
on 17 Jun 2013

please read listing of the function ppval (in Command Window: >> open ppval).

b = pp.breaks;
c = pp.coefs;
k = pp.order;
xx = linspace(b(2),b(3),200);
sizexx = size(xx);
lx = numel(xx); 
xs = reshape(xx,1,lx);
[~,index] = histc(xs,[-inf,b(2:l),inf]);
xs = xs-b(index);
v = c(index,1);
for i=2:k
   v = xs(:).*v + c(index,i);


v = ppval(pp,xx);

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply today