4
$\begingroup$

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?

$\endgroup$
6
  • 5
    $\begingroup$ Providing some code for f[t] or using Integrate might help :) $\endgroup$
    – Öskå
    Aug 7, 2014 at 9:55
  • 3
    $\begingroup$ @Roman The question seems closer to mathematics than related to Mathematica (I am also curious about this I just think It will have a better chance at Mathematics stack exchange). $\endgroup$ May 1, 2023 at 15:04
  • $\begingroup$ Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Mathematica Meta, or in Mathematica Chat. Comments continuing discussion may be removed. $\endgroup$
    – Kuba
    May 2, 2023 at 5:06
  • $\begingroup$ What you are asking is impossible in general. You only know the value of the function at $t_1,...,t_n$ and yet the integral is over $[0,\infty)$. Unless you know something about the values of the function for $t>t_n$ you can't find the integral over the infinite interval $[t_n,\infty)$ and hence over $[0,\infty)$. Interpolation is relatively easy, while extrapolation is very hard. $\endgroup$
    – Somos
    May 2, 2023 at 23:43
  • $\begingroup$ Related: mathematica.stackexchange.com/questions/162768/… -- though I find the time step being slightly greater than twice the period of the carrier frequency a bit troubling.... $\endgroup$
    – Michael E2
    May 3, 2023 at 3:49

3 Answers 3

4
$\begingroup$

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).

$\endgroup$
4
+500
$\begingroup$

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

among maybe many others.

$\endgroup$
2
  • $\begingroup$ Nice! Your writeup clarifies a few ideas and is exactly the kind of discussion I wanted to stimulate with the bounty. The finite summation error is, in my opinion, not worthy of much more discussion, as it is covered by most discussions on the techniques of numerical integration. The truncation error is really where the interest lies. Great work, thanks! $\endgroup$
    – Roman
    May 7, 2023 at 17:11
  • $\begingroup$ @Roman I placed the truncation error section first instead of last and reduced the summation error to a list of possible methods. It does seem to make it easier to read $\endgroup$ May 8, 2023 at 5:15
3
$\begingroup$

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]]}.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.