# High precision numerical integration of discrete data with Matlab

I have discrete data of a function plotted below: The "Y" values of the function near "X=1.57" are very close to each other and zero, like 9.25558265263186E-11 and 5.92357284527227E-11.

Using "Trapz" function directly, or using "Integral" function indirectly as

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


fails with more than 30% relative error. I could have more data points for the function, however, using 100 points or 1 million points does not change the result.

My data for the function and the reference solution are from a physical problem calculated through multiple long computational steps(very long to explain). Since my code and the the computations have already been verified for a *better function, and since the same problem was reported in a paper (that singularity begins for this setting and the numerical integration becomes harder), I think the problem is related to the numerical integration, probably because of the roundoff errors and the 16digits precision of Matlab.

I also did all of the computations from the first of the simulation using VPA. for each line of code, i used VPA atleast once (In Matlab website, If i understand correctly, it is written that if one of the numbers/variables are VPA, whole computation in that line is done using higher precision). I mean, my data (y values) were calculated using VPA before being given to TRAPZ function. Then, I used

trapz(x,vpa(y))% y also has already been calculated using VPA.


About X values, i could not make it higher precision I am working on that.

With and without VPA: 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.

I wrote the same question in another website, and people used my data (*for 100k) points and also got the same result (0.29...). They explained me that such 30% error is not happening because of the integration and it is related to my data and reference solution. Is there any chance that we missed some point here?

• in my physical problem, I can change some parameters, which leads to a better function (not near-zero values around x=1.57). The more my function becomes like this, the more my error increases.

• 100k_points.csv

• Welcome to Scicomp. Are you looking for a bug in your current implementation, or a better numerical integration method for your problem? Could you give us the Function you plottet in mathematical form? That way someone might be able to reproduce the problem for themselves. Dec 31, 2021 at 13:47
• @MPIchael thanks. Actually maybe someone knows if that 30% relative error for this plotted function should be expected? Is this function bad enough to expect that error? I don't have any mathematical function. It is evaluated through multiple long computational steps. Since It has been verified for other settings/parameters, and as the same problem was reported by a paper, I guess my problem might be related to integration.
– tio
Dec 31, 2021 at 17:04
• Why are you interpolating the function and using trapz to integrate- why not simply sum the data values times $\Delta x$? Dec 31, 2021 at 17:56
• Trapz does the same similar to the sum. I used the interpolation as a second approach.
– tio
Dec 31, 2021 at 18:00

Since it looks like your discretized function is monotonically decreasing, you can calculate an effective upper and lower bound of the numerical integration using left and right Riemann sums.

There is a way to compute accurate summations of floats which avoids repeated summation truncation error. I don't know if Matlab has an implementation of this function or not, but in Python the built-in math.fsum function implements this. Another alternative implementation is to use Kahan summation.

Using the accurate floating point summation, I computed the left and right Riemann sums using:

import numpy
import math
x = data[:,0]
y = data[:,1]
dx = x[1:] - x[:-1]
left_area = y[:-1] * dx
right_area = y[1:] * dx
left_int = math.fsum(left_area)
right_int = math.fsum(right_area)


This gives the bounds for numerical integration error as:

min value: 0.2996090880138866
max value: 0.29962515632195036
delta: 1.6068308063887926e-05


From these results along with the result produced by trapz or even simps, I think it's safe to say that the numerical integration is accurate to at least 4 digits, so the problem is either in your simulation computing the curve or how you calculated the expected integral value.