Finding the area under a polynomial equation in MATLAB

Asked by mr mo

mr mo (view profile)

on 30 Apr 2018
Latest activity Commented on by mr mo

mr mo (view profile)

on 30 Apr 2018
Accepted Answer by John D'Errico

John D'Errico (view profile)

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

John D'Errico (view profile)

on 30 Apr 2018
Edited by John D'Errico

John D'Errico (view profile)

on 30 Apr 2018

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

mr mo (view profile)

on 30 Apr 2018
Thanks a lot. This works very well.

Answer by James Tursa

James Tursa (view profile)

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

mr mo (view profile)

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

James Tursa (view profile)

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

mr mo (view profile)

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