4
$\begingroup$

Only the first two sections of this long question are essential. The others are just for illustration.

Background

Advanced quadratures such as higher-degree composite Newton–Cotes, Gauß–Legendre, and Romberg seem to be mainly intended for cases where one can finely sample the function but not integrate analytically. However, for functions with structures finer than the sampling interval (see Appendix A for an example) or measurement noise, they cannot compete with simple approaches such as the midpoint or trapezoid rule (see Appendix B for a demonstration).

This is somewhat intuitive as, e.g., the composite Simpson rule essentially “discards” one quarter of the information by assigning it a lower weight. The only reason such quadratures are better for sufficiently simple functions is that properly handling border effects outweighs the effect of discarded information. From another point of view, it is intuitively clear to me that for functions with a fine structure or noise, samples that are remote from the borders of the integration domain must be almost equidistant and have almost the same weight (for a high number of samples). On the other hand, quadrature of such functions may benefit from a better handling of border effects (than for the midpoint method).

Question

Assume that I wish to numerically integrate noisy or fine-structured one-dimensional data.

The number of sampling points is fixed (due to function evaluation being costly), but I can freely place them. However, I cannot place sampling points interactively, i.e., based on results from other sampling points. I also do not know potential problem regions beforehand. So, something like Gauß–Legendre (non-equidistant sampling points) is okay; adaptive quadrature is not since it requires interactively placed sampling points.

  • Have any methods going beyond the midpoint method been suggested for such a case?

  • Or: Is there any proof that the midpoint method is best under such conditions?

  • More generally: Is there any existing work on this problem?

Appendix A: Specific example of a fine-structured function

I wish to estimate $\int_0^1f(t)\, \mathrm{d}t$ for: $$ f(t) = \sum_{i=1}^{k} \frac{\sin(ω_i t-φ_i)}{ω_i},$$ with $φ_i∈ [0,2π]$ and $\log{ω_i} ∈ [1,1000]$. A typical function looks like this:

superimposed sines

I chose this function for the following properties:

  • It can be integrated analytically for a control result.
  • It has fine structure on a level that makes it impossible to capture all of it with the number of samples I am using ($<10^2$).
  • It is not dominated by its fine structure.

Appendix B: Benchmark

For completeness, here is a benchmark in Python:

import numpy as np
from numpy.random import uniform
from scipy.integrate import simps, trapz, romb, fixed_quad

begin = 0
end   = 1

def generate_f(k,low_freq,high_freq):
    ω = 2**uniform(np.log2(low_freq),np.log2(high_freq),k)
    φ = uniform(0,2*np.pi,k)
    g = lambda t,ω,φ: np.sin(ω*t-φ)/ω
    G = lambda t,ω,φ: np.cos(ω*t-φ)/ω**2
    f = lambda t: sum( g(t,ω[i],φ[i]) for i in range(k) )
    control = sum( G(begin,ω[i],φ[i])-G(end,ω[i],φ[i]) for i in range(k) )
    return control,f

def midpoint(f,n):
    midpoints = np.linspace(begin,end,2*n+1)[1::2]
    assert len(midpoints)==n
    return np.mean(f(midpoints))*(n-1)

def evaluate(n,control,f):
    """
    returns the relative errors when integrating f with n evaluations
    for several numerical integration methods.
    """
    times = np.linspace(begin,end,n)
    values = f(times)
    results = [
            midpoint(f,n),
            trapz(values),
            simps(values),
            romb (values),
            fixed_quad(f,begin,end,n=n)[0]*(n-1),
        ]

    return [
            abs((result/(n-1)-control)/control)
            for result in results
        ]

method_names = ["midpoint","trapezoid","Simpson","Romberg","Gauß–Legendre"]

def med(data):
    medians = np.median(np.vstack(data),axis=0)
    for median,name in zip(medians,method_names):
        print(f"{median:.3e}   {name}")

print("superimposed sines")
med(evaluate(33,*generate_f(10,1,1000)) for _ in range(100000))

print("superimposed low-frequency sines (control)")
med(evaluate(33,*generate_f(10,0.5,1.5)) for _ in range(100000))

(I here use the median to reduce the influence of outliers due to functions that have only high-frequency content. For the mean, the results are similar.)

The medians of the relative integration errors are:

superimposed sines
6.301e-04   midpoint
8.984e-04   trapezoid
1.158e-03   Simpson
1.537e-03   Romberg
1.862e-03   Gauß–Legendre

superimposed low-frequency sines (control)
2.790e-05   midpoint
5.933e-05   trapezoid
5.107e-09   Simpson
3.573e-16   Romberg
3.659e-16   Gauß–Legendre

I already asked this question two months ago on Computational Science SE and placed a bounty on it without success so far.

$\endgroup$
  • $\begingroup$ So...[obligatory disclaimer that I don't know anything about anything] would it help to have a way to "cumulatively" perform the higher order newton-cotes methods? By which I mean, given the starting point and the signal, performing the integral based upon past points each "time step", that way the integrating function just uses the number of points it has? (I'm sorry, I am just having trouble understanding the question.) $\endgroup$ – heather Jul 12 '18 at 1:23
  • $\begingroup$ @heather: I do not understand your suggestion and can only guess that what you are suggesting is a composite Newton–Cotes method. In that case: no it doesn’t help (trapezoid and Simpson in the above are the composite Newton–Cotes methods of order 2 and 3. Can you elaborate your troubles understanding the question, so I can try to address them? $\endgroup$ – Wrzlprmft Jul 12 '18 at 5:46
  • $\begingroup$ no, I'm not quite talking about composite Newton-Cotes methods...I'm talking about versions of these that do not require knowledge of the function and allow for time step variance. I really apologize; I'm not very good at conveying what I mean. Feel free to ping me in some chatroom to talk. $\endgroup$ – heather Jul 16 '18 at 17:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.