Asked by mr mo
on 30 Apr 2018

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 by John D'Errico
on 30 Apr 2018

Edited by John D'Errico
on 30 Apr 2018

Accepted Answer

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.

Answer by James Tursa
on 30 Apr 2018

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
on 30 Apr 2018

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
on 30 Apr 2018

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
on 30 Apr 2018

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.