vpaintegral of discrete data

I have discrete data of the attached function stored in the array 'x' and 'y'.

I can have infinite number of points for this function.

First, I used the Trapezoidal rule.

Second, I used the interpolation function of Matlab to be able to use the built-in Integral function of Matlab as follows:

f1 = griddedInterpolant(x,y); %I used also spline and other methods
f2 = @(t) f1(t);
Result = Integral(f2,0,pi)

Comparing to my reference solution, both of those methods provide significant error. (I verified my code and the way i achieved the reference solution with another smooth function). The difference between using 1000 points or 1million points is negligible. I guess this shows that the problem is not related to round off error. However, I wanted to also use VPAINTEGRAL, but i cannot do it because i cannot convert the interpolated function to a symbolic form. Appreciate any suggestion how to use VPAINTEGRAL (If it can help...)

[Torsten]

What in your opinion is a "significant error" ?

[Msimon q]

This function was numerically integrated by someone in a paper with significantly lower error. I am sure that something is wrong.

Answer 1

[John D'Errico]

Sorry. vpaintegral CANNOT be used. It simply does not apply to data points, as you then have no function. And the high, multi-digit precision of vpa is then meaningless.

Could you form a spline approximation to the data sequence, and then integrate that using vpaintegral? Well, technically yes, but that would be no better than using integral.

To go back to an early comment of yours, there is a very good chance that you used interpolation wrong, when you tried to then use integral. If so, then it was no better than using trapz.

I'll show an example, but since you posted only a picture of the function instead of the function itslef, I cannt really help you. At best, I can show an example. I'll integrate the sine functino, from 0 to pi. The theoretical value is 2.

syms X
int(sin(X),[0,pi])

ans = 

2

Can we form a good approximation using trapz? Note that here, trapz will tend to under-estimate the true value. I'll not take that large of a sample, so we can see the difference between the different approaches I take.

fun = @sin;
x = linspace(0,pi,20);
y = fun(x);

format long g
trapz(x,y)

ans = 
          1.99544131832019

Can we now use integral? Again, BE CAREFUL! For example, if I do it in a very simplistic way, without taking care...

interpfun = @(u) interp1(x,y,u);
integral(interpfun,0,pi,'abstol',1e-14)

ans = 
          1.99544121754747

This was little different from trapz. However, try this:

splinterpfun = @(u) interp1(x,y,u,'spline');
integral(splinterpfun,0,pi,'abstol',1e-14)

ans = 
          1.99999851051087

So when we do the interpolation correctly to be more accurate, now the integral is much more close to the true value. Not perfect. But then we cannot easily attain perfection, if all we are willing to provide are data points.

Again, vpaintegral is essentially useless here. (Could you use it? SIGH. With some effort and a far deeper understanding of splines than I am willing to teach here, you could, but it would do no better than using integral itself as I did.)

Looking back at your question, I do see you used griddedInterpolant. As it turns out, the difference is zero. Interp1 and griddedInterpolant are the same here.

[Msimon q]

Thanks for your explanation. But I don't get "there is a very good chance that you used interpolation wrong.." because i wrote what i did. It is pretty same as what you wrote. I also used 'abstol',1e-16 but i did not help...

my "y" array has some values like this:

0.267955045816377 0.267955045817655 0.267955045818781

Since they are very close to each other, i guess that both Trapz and Integral function do some kind of roundoff error.

[Torsten]

Let's stop guessing.

Please show your x- and y-data, your code, the integral you get with your code and the integral you expect to get.

[Msimon q]

Thanks. For simplicity and because of the symmetry, I use first half of the curve for another similar curve. I attach the curve data for two cases:10 points and 100k points. The first column is the X axis.

Using 10 points and Trapz, I get : 0.299517049200365 , (I get nearly the same with "Integral" function as i explained above.)

Using 100k points and Trapz, I get: 0.299617122167918,

The reference solution is : 0.220712757091533,

For a better curve (where the values of "Y" around "X=1.5" is not like this one), my numerical integration matches the corresponding reference solution perfectly.

Thanks for your help

[Torsten]

Such a big difference in the value of the integral cannot be explained by the integration method or the number of points.

There must be a substantial difference in the functional equation for your data and for the reference data.

[Jan]

See your other question: The only reasonable explanation is, that the reference solution "0.2207" does not match the shown data. Where does this reference come from?

[Msimon q]

@Torsten, Thanks https://de.mathworks.com/matlabcentral/answers/1619300-numerical-integration-of-discrete-data-using-higher-precision#comment_1910655