How to integrate a function which is only known at discrete points?

Question

I have an integration to do. I want to integrate

$$ \int_0^\infty \sin^2(2\pi t)f(t)\mathrm{d}t $$

where $f(t)$ is known only at discrete values given in an array in the form $\{t_i,f_i\}$ with $i=1\ldots n$.

The time steps in the array is 1.1s. Can you please suggest a method to do this? I tried using the Trapezoidal method for numerical integration but gave a bad approximation. Is there an easy method with inbuilt function or another method?

Answer 1

If you were to have, for example,

dt = .01;
tbl = Table[{t, Exp[Cos[t]]}, {t, 0, 10, dt}];

(that is, your $f(t)$ corresponds to my tbl) then another way is

Total @ MapThread[
    Sin[2*Pi #1]^2 * #2 &,
    Thread @ tbl
   ] * dt

But of course any technique can be easily used (trapezoidal or more sophisticated approaches).

Answer 2

Consider the following example:

integrand[x_] := Sin[2*x]/(1 + x^2);

Computing the integral of that function from 0 to infinity directly has 2 sources of errors:

• Truncation error by taking a $x_{\textit{max}}$: $\int_0^\infty \mathrm{integrand}(x) \mathrm{d} x \longrightarrow \int_0^{x_{\mathrm{max}}} \mathrm{integrand}(x) \mathrm{d} x$
• finite summation error: $\int_0^{x_{\mathrm{max}}} \mathrm{integrand}(x) \mathrm{d} x \longrightarrow \sum ...$

In the following, I will focus on the truncation error. At the end of this post, I will provide a list of resources for the finite summation error.

Outline:

• Improving the truncation error with summation accelerators

  • The method
  • The code

• Methods to minimize the finite summation error

Improving the truncation error with summation accelerators

I do not know what the name of this is but there is a method that consists of replacing an integral of an oscillating integrand with a heavy tail by an alternating series and then using a series accelerator. I have a vague memory that the usual series accelerators do not help much when the series already converges quickly but I might be wrong. In any case, for the integrand considered here, I got a better result (that said I did not test this extensively for different parameters to dismiss the scenario of a fluke).

The method

Let $S$ be the set of all zeroes of the integrand. If we consider that the integrand always changes sign at its zeros in the domain of integration, then we may re-write the integral as:

$\int_0^\infty \mathrm{integrand}(x) \mathrm{d} x= \sum_{z_i \in S} \int_{z_{i}}^{z_{i+1}} \mathrm{integrand}(x) \mathrm{d} x=\sum_i (-1)^i h_i$

The alternating sign is due to the fact that the function has a fixed sign between its zeros and changes sign at each crossing.

We can then apply the Shanks accelerator to the partial sums.

This is what is done below.

Code

Note: •=[Bullet] and ⎵=[UnderBracket]

Simpson integration:

•SimpsonIntegration[pts_,step⎵size_]:=
Module[{weights},
weights=
Join[{1}
    ,
    ArrayPad[{4,2}
            ,
            {0,Length[pts]-4}
            ,
            "Periodic"
    ]
    ,
    {1}
]
;
pts.weights*(step⎵size/3)
]

Data points:

max = 100; n = max*100; min = 0; data = 
 Subdivide[min, max, n] // Map[{#, integrand[#]} &];

The function that multiplies the $\sin$ function does not have any zeroes in the domain of integration. Hence, the zeroes are those of $\sin 2x$ at $k \pi/2 $. Thus, points should be grouped into intervals $[\frac{k\pi}{2},\frac{(k+1)\pi}{2}]$

parts = GatherBy[data, Floor[2*#[[1]]/Pi] &];

In the data points, the boundary points of each interval, that is the points ${x=k\pi/2,\mathrm{integrand}(x)=0}$, are not necessarily included. Adding these points seems to lead to better accuracy which is what is done below. Note that adding these points will break the uniformity of the distribution of points in general as the distance from neighbors can be smaller than the original step size. However, the formula for the Simpson rule for non-uniform step sizes at https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data shows that the non-uniform step size contributions would be multiplied by the value of the integrand at the zeros, that is, the contributions from non-uniformity are all zero in this case.

c = 0;
parts•bulk = 
  parts[[2 ;; -2]] // DeleteCases[{_, 0}] // 
   Map[(c++; Join[{{c*Pi/2, 0}}, #, {{(c + 1)*(Pi/2), 0}}]) &];
c++; parts = 
 Join[{Join[parts[[1]] // DeleteCases[{Pi/2, _}], {{Pi/2, 0}}]}, 
  parts•bulk, {Join[{{c*Pi/2, 0}}, 
    parts[[-1]] // DeleteCases[{c*Pi/2, _}]]}]; 

Now we calculate the integrals at each subdivision:

[Edit written months after: When using Simpson's rule, the number of points should be odd which I believe was not the case here but the result obtained is still fairly close to the reference]

integrals = 
  parts // 
   Map[•SimpsonIntegration[#[[All, 
        2]], #[[3, 1]] - #[[2, 1]]] &]; 

The direct sum without the series acceleration:

integrals // Total // N

(* 0.515875 *)

The reference:

Integrate[integrand[x], {x, 0, Infinity}] // N

(* 0.515905663339148` *)

Now we consider accelerating the series. First, we compute partial sums:

partial⎵sums = N@Accumulate@integrals;

Then we define the Shanks formula from here:

Shanks[A_, 
   n_] := (A[n + 2] A[n] - A[n + 1]^2)/(A[n + 2] + A[n] - 2 A[n + 1]);

and we use this on the last few partial sums:

Shanks[partial⎵sums[[#]] &, 
 Length[partial⎵sums] - 2]

(* 0.5159066435204446`  *)

Methods to minimize the finite summation error

• Low-order composite Newton-Cotes

• rational approximations (for example https://www.sciencedirect.com/science/article/pii/S0377042713001714)

• sinc interpolation(for example https://arxiv.org/abs/1902.08562) (maybe apply this only to the function that multiplies the $\sin$ in the integral.)

• Polynomial least squared fit? Did not try but see https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas

among maybe many others.

Answer 3

Let's denote values {t, f(t)} as F, then interpolate this array with ListInterpolation.

fx = ListInterpolation[F[[All, 2]], {F[[1, 1]], F[[-1, 1]]}]

Now we may use Integrate with fx

Integrate[Sin[2*Pi*t]*fx[t], {x, F[[1, 1]], F[[-1, 1]]}]

As you see, it's not exactly what you want: the domain of integration is {F[[1, 1]], F[[-1, 1]]}.