Numeric integration - singularity

Question

I have discrete data points of a function(I can have infinite number of points of the function), which is something like this:

I used trapezoidal function of Matlab for a numerical integration.

The problem is that increasing the number of points, does not reduce the error of the integration. As someone already answered, the problem lies on the spikes. But,

• why increasing the number of points does not decrease the error? Using 10k points or 10million points makes no difference.

Isn't really 10 million points enough to capture all details of this function?

Answer 1

I have had this problem with MATLAB before too, the problem lies in how it stores numbers. Theoretically the error would decrease but actually you are summing up more regions, each of which has a relatively increasing error. Here's an example:

---

Say I want to work out the derivative numerically at a given point, we have: $$\frac{f(x+h)-f(x)}{h}$$ however, when making $h$ very small it will be fine if its a simple number (like $10^a)$ but when it has a large number of non-zero digits, matlab will store it rounded after the output, and so the value of $f(x+h)$ will not have the precision you want. Also, during the summation process there will be rounding after every addition.

---

Matlab by default uses 16 digits, however if you use the vpa function you can increase it, so you could write your own code for numerical integration that incorporates this, the number of digits vpa used can also be changed.

vpa Variable precision arithmetic. R = vpa(S) numerically evaluates each element of the double matrix S using variable precision floating point arithmetic with D decimal digit accuracy, where D is the current setting of DIGITS. The resulting R is a SYM.

vpa(S,D) uses D digits, instead of the current setting of DIGITS.
D is an integer or the SYM representation of a number.

It is important to avoid the evaluation of an expression using double
precision floating point arithmetic before it is passed to vpa.
For example,
   phi = vpa((1+sqrt(5))/2)
first computes a 16-digit approximation to the golden ratio, then
converts that approximation to one with d digits, where d is the current
setting of DIGITS.  To get full precision, use symbolic arguments,
   phi = vpa((1+sqrt(sym(5)))/2)
or
   s = sqrt(sym(5))
   phi = vpa((1+s)/2)

Additional examples:
   vpa(pi,780) shows six consecutive 9's near digit 770 in the
      decimal expansion of pi.

   vpa(hilb(2),5) returns

      [    1., .50000]
      [.50000, .33333]

See also double, digits, subs.

Documentation for vpa
Other functions named vpa

You can recall this yourself using the command

help vpa