Finding the area under a polynomial equation in MATLAB

Hi. I have two vectors like these.

 X=[0   0.0300   0.0333   0.0469   0.0575   0.0649   0.0960   0.1262   0.1610   0.1717   0.2017   0.2411   0.2711   0.3000   0.3035   0.3368   0.3711   0.4068   0.4411   0.4758   0.5253   0.5553   0.5853   0.6000]
 Y=[0 126.2766 140.0068  190.5906  211.6557  219.3601  236.8709  253.2271  266.4663  269.3085  273.6913  279.3095  283.5479  286.1191  286.2725  282.7424  277.6845  272.1010  266.4556  260.1458  248.5416  241.1381  233.2857  230.3781]

I use the polyfit command to fit a polynomial curve of degree 6 to my data-set.

 B=polyfit(X,Y,6)

Now my question is how can I calculate the area under this polynomial curve?

Thanks a lot.

Answer 1

[John D'Errico]

How do you doit? Basic calculus?

X=[0   0.0300   0.0333   0.0469   0.0575   0.0649   0.0960   0.1262   0.1610   0.1717   0.2017   0.2411   0.2711   0.3000   0.3035   0.3368   0.3711   0.4068   0.4411   0.4758   0.5253   0.5553   0.5853   0.6000];
Y=[0 126.2766 140.0068  190.5906  211.6557  219.3601  236.8709  253.2271  266.4663  269.3085  273.6913  279.3095  283.5479  286.1191  286.2725  282.7424  277.6845  272.1010  266.4556  260.1458  248.5416  241.1381  233.2857  230.3781];

xinterp = linspace(min(X),max(X),100);
P = polyfit(X,Y,6);
plot(X,Y,'ro',xinterp,polyval(P,xinterp),'b-')

As models go, I'd call it your basic piece of crapola. And DON'T use a higher order polynomial, thinking it will do better. Just use a spline instead, if you ABSOLUTELY insist on a model. Actually, pchip (shape preserving interpolant) would be my choice here.

If I had to guess, that looks like essentially a piecewise linear curve, with 5 linear segments. But why do you need to integrate a poorly fit polynomial at all?

The area under the curve is...

trapz(X,Y)
ans =
     149.16

In fact, this is arguably a better estimator of the area than the integral of that crappy polynomial model.

What do I get for the pchip integral? My SLM toolbox has a nice little utility to integrate a spline. But it is easily enough done without that tool.

pp = pchip(X,Y);
integral(@(x) ppval(pp,x),min(X),max(X))
ans =
       149.28

That is not too far off from what trapz gave. Bother are about as good as you can get.

Can we integrate the polynomial model? Again, of course. As I said, basic calc suffices. (Or, you could use polyint.)

P = polyfit(X,Y,6);
Pint = [P./(7:-1:1),0];
diff(polyval(Pint,[min(X),max(X)]))
ans =
     149.78

Or if you are feeling lazy,

integral(@(x) polyval(P,x),min(X),max(X))
ans =
     149.78

In any event, I'd use what pchip or what trapz gave here as better estimators.

[mr mo]

Thanks a lot. This works very well.

Answer 2

[James Tursa]

Just integrate the polynomial and evaluate at the endpoints. E.g.,

IB = polyint(B);
area_under_poly = polyval(IB,X(end)) - polyval(IB,X(1));

[mr mo]

Thanks but if I have these vectors, the calculated area must be 102.5 but this code return 112.5

 x=[0    10    15]
 y=[0    10    11]
 B=polyfit(x,y,2)

[James Tursa]

By hand (using the fact that x(1)=0):

>> (1/3)*B(1)*x(end)^3 + (1/2)*B(2)*x(end)^2 + B(3)*x(end)
ans =
  112.5000

Why do you think the answer must be 102.5?

[mr mo]

Because it is obvious and can be seen by using the plot command

 plot(x,y,'-rs')
 grid on

if the vectors are like these

 x=[0    10    15]
 y=[0    10    15]

then the area will be equal to 112.5