High precision numerical integration of discrete data with Matlab

Question

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

Answer 1

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
data = numpy.loadtxt('100kPoints.csv', delimiter=',')
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.