about trapz and quad integral

Hi all;

I have the folowing data:

x_1=data(:,1);

y_1=data(:,2);

xx=linspace(x(1),x(end),10000);

pp=spline(x,y);

yy=ppval(pp,xx);

%%integral over [x(1),x(end)]

I=trapz(xx,yy)

Int1=quad(@ppval,x_1(1),x_1(end),[],[],pp)

the value of the integral calculated using quad and trapz are different!!!!!!

I=trapz(x,y) gives I=1.353085426296705e+09

and Int1=2.329171651204020e+10 and I obtained the following warning:

Warning: Maximum function count exceeded; singularity likely.

Why there is such differnce ?

how can I avoid this warning?

thank yoy for your help

cheers

[Yu Jiang]

Hi Adam

I think the question can be better answered if the following issues can be explained:

• First, the code uses trapz (See Documentation) to integrate over the interval of [x(1),x(end)] but use quad (See Documentation) over the interval of [x_1(1),x_1(end)]. Can you confirm that x and x_1 are the same?
• Can you share your data set, i.e., the matrix data? Or can you let me know how I can generate this data set to reproduce the issue?

- Yu

[Adam]

Hi Yu,

thanks for your reply.

Find attached please the data: data=psf110

cheers

Answer 1

[Yu Jiang]

Let us first take a look at the differences between quad (See Documentation) and trapz (See Documentation). quad approximates the integral of a given function by using the quadrature method, i.e., breaking the area down into rectangles. On the other hand, trapz approximates integral by breaking the area down into trapezoids. Therefore, trapz generally provides more accurate results.

From the data set you provided, a huge spike around x=0 can be observed. This seems to be related with the warning message “singularity likely.” According to the diagnostic section in the documentation , 'Maximum function count exceeded' indicates that the integrand has been evaluated more than 10,000 times. A nonintegrable singularity is likely. There seems to be no way to change that maximum number.

If you would like to calculate this integral by using quadrature method, there are two workarounds as follows.

• Note that instead of using trapz, which uses recursive adaptive Simpson quadrature, you may want to use the function integral (See Documentation) , which uses global adaptive quadrature. In this case, the last line in your code can be changed to

Int1 = integral(@(x) ppval(pp,x),xx(1),xx(end))

I have tested it. The result is 1.3307e+09.

• If quad has to be used and you would like to improve the accuracy and avoid the warning, I would suggest you break down the problem into multiple integration problems over smaller intervals, calculate them individually, and sum up the results.