11
$\begingroup$

I mean $$\int_ 0^{2\pi}\frac {\sin (10050 t)\sin (10251 t)\cos (2022 t)} {\sin (50 t)\sin (51 t)} dt.$$

Here is my unsuccessful trial:

NIntegrate[Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t], {t, 0,2*Pi}, AccuracyGoal -> 5, PrecisionGoal -> 5]

NIntegrate::ncvb : NIntegrate failed to converge to prescribed accuracy after 9
recursive bisections in t near {t} = {3.15386} . NIntegrate obtained 11.171230707439639 and 8.934855726568145 for the integral and error estimates .

I don't find an answer in the documentation (see the "Finite Region Oscillatory Integration" section).

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Try TrigReduce on the integrand $\endgroup$
    – Andreas
    Commented Sep 27, 2022 at 8:26
  • 1
    $\begingroup$ I am interested rather in its numeric calculation (The documentation presents examples of such type.). $\endgroup$
    – user64494
    Commented Sep 27, 2022 at 10:54
  • 3
    $\begingroup$ Maybe MinRecursion -> 8, Method -> {"GaussKronrodRule", "Points" -> 21}. $\endgroup$
    – Goofy
    Commented Sep 27, 2022 at 11:30
  • 1
    $\begingroup$ @Goofy: Many thanks from me to you. This is it. Can you present your comment as an answer, elaborating it? In particular, how did you come to these options? $\endgroup$
    – user64494
    Commented Sep 27, 2022 at 11:51
  • 1
    $\begingroup$ It should be noticed that NIntegrate[ Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t], {t, 0, 2*Pi}, Method -> {"LevinRule", "Kernel" -> Cos[2022*t]}, MaxRecursion -> 8] results in 8.75152 which is far away from the right one. $\endgroup$
    – user64494
    Commented Sep 27, 2022 at 13:01

6 Answers 6

Reset to default
7
$\begingroup$

Here's a fast way:

NIntegrate[
  Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t], {t, 0, 
   2*Pi}, MinRecursion -> 8, 
  Method -> {"GaussKronrodRule", "Points" -> 21}] // AbsoluteTiming

(*{0.224591, 18.849555920879993`}*)

The integrand is analytic and the Gauss rule converges fast on analytic functions, so a high-order rule is probably going to be helpful. There are about 44K zeros. The Gauss-Kronrod rule is order 3*21+2 and the interval is divided into 2^8 subintervals. Together that's 2^8*65 or almost 17K zeros. Since the result is 5.999999999790308 * Pi, something must be working well. It seems like undersampling, but aliasing is a possible explanation for the unexpected accuracy (which I cannot explain).

SolveValues[
  Cos[2022*t]*Sin[10050*t]*Sin[10251*t] == 0 && 0 <= t <= 2*Pi, 
  t] // Length

(*44245*)

NIntegrate[Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t],
  Evaluate@Flatten@{t, 
     SolveValues[
      Cos[2022*t]*Sin[10050*t]*Sin[10251*t] == 0 && 0 <= t <= 2*Pi, 
      t]}, MinRecursion -> 0, 
  Method -> "GaussKronrodRule"] // AbsoluteTiming

(*{19.1047, 18.849555921538766`}*)

The result is 6.000000000000003 * Pi, which is more accurate but takes much longer.

Improve this answer
$\endgroup$
7
  • $\begingroup$ +1. I accept your answer, but others also are good and deep. $\endgroup$
    – user64494
    Commented Sep 29, 2022 at 5:28
  • 1
    $\begingroup$ It would be interesting to know if this would be as effective if the big peak was at unknown location, rather than at the start and end. $\endgroup$
    – mikado
    Commented Sep 29, 2022 at 16:30
  • $\begingroup$ Perhaps you saw my answer where I used Reduce to find the zeroes. I did not check if some of the other solutions from the output of Reduce that had auxiliary constants were distinct zeroes from the ones I used. I corrected my answer and checked that we found the same zeroes after using sort. Reduce is about 7 times faster for finding the zeroes but tedious to use to obtain numerical values. $\endgroup$
    – userrandrand
    Commented Sep 29, 2022 at 17:26
  • 1
    $\begingroup$ @mikado there is a peak in the middle but Plot misses it if the number of points is not increased from the default choice. For example Plot[Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t], {t, 0, 2*Pi}, PlotRange -> All, PlotPoints -> 4000]. $\endgroup$
    – userrandrand
    Commented Sep 29, 2022 at 18:16
  • 1
    $\begingroup$ @mikado changing the upper bound to RandomReal[{1, 4}]*Pi/2 , MinRecursion to 10 and MaxRecursion to 20 (for example), using Do, after evaluating the integral 100 times, I found an average of 0.1 seconds per integral with no error using the method above. $\endgroup$
    – userrandrand
    Commented Sep 29, 2022 at 18:26
7
$\begingroup$

Note: As OP has pointed out in a comment, a similar calculation is available here on Math SE.

OPs integral can be written as $$ \int_0^{2\pi} f(50t) f(51t) \cos(2022 t) dt $$ with the auxiliary function

f[x_] := Sin[201*x]/Sin[x];

The division can be carried out explicitly, because (X^201-Y^201)/(X-Y) is a polynomial. The answer is $$ f(x) = e^{-200 ix} + e^{-198 ix} + \ldots + e^{-2ix} + 1 + e^{2ix} + \ldots + e^{198 ix} + e^{200 ix} $$ or as Mathematica code:

falternative[x_]:=Sum[Exp[I*k*x],{k,-200,200,2}];

To check this, use

f[x]-falternative[x]//FullSimplify
(* 0 *)

OPs integrand is therefore

integrand = falternative[50*t]*falternative[51*t]*Cos[2022*t]//TrigToExp//Expand;

This returns a linear combination of terms of the form $e^{i \omega t}$ with integer $\omega$, and the coefficient of $\omega = 0$ is $3$, therefore the value of the integral is $$ \int_0^{2\pi} f(50t) f(51t) \cos(2022 t) dt = 6\pi $$

Improve this answer
$\endgroup$
2
  • 4
    $\begingroup$ +1. Thank you for your work. I knew the exact answer from https : // math . stacketchange . com/questions/4508185/evaluate - int2 - pi - 0 - frac - cos2022t - sin10050t - sin10251t - sin50t - sin #4508185 and that was a reason to ask it. I am interested rather in its numeric calculation (The documentation presents examples of such type.). $\endgroup$
    – user64494
    Commented Sep 27, 2022 at 10:53
  • 1
    $\begingroup$ Really interesting question and answer. EE's sometimes call sin(x)/x the "sync" functions, and there are variants based on pi, etc. $\endgroup$
    – JosephDoggie
    Commented Sep 27, 2022 at 19:26
4
$\begingroup$

With patience, Monte Carlo methods can be effective

NIntegrate[
 Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t], {t, 0, 
  2*Pi}, AccuracyGoal -> 5, PrecisionGoal -> 5, MaxRecursion -> 10^6, 
 MaxPoints -> 10^8, Method -> "AdaptiveQuasiMonteCarlo"]
(* 18.8496 *)

% - 6 π
(* 2.4873*10^-8 *)
Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ This is good, but there is no need for Monte-Carlo -- just use "GlobalAdaptive" without symbolic pre-processing and without singularity handling. $\endgroup$
    – Anton Antonov
    Commented Sep 28, 2022 at 23:55
  • $\begingroup$ +1. Unexpected accuracy of Monte-Carlo method. $\endgroup$
    – user64494
    Commented Sep 29, 2022 at 5:30
  • $\begingroup$ @AntonAntonov there are reasons not to use Monte Carlo as the last method you try for a difficult integral, but it is often a very good first choice. $\endgroup$
    – mikado
    Commented Sep 29, 2022 at 18:25
  • $\begingroup$ I have never heard that rule of thumb for 1D integrals. :) (Only for dimensions larger than, say, 5...) $\endgroup$
    – Anton Antonov
    Commented Sep 29, 2022 at 18:34
  • 2
    $\begingroup$ @AntonAntonov in principle, I am interested in those things. A major interest is in estimating the effectiveness of electronic systems (and optimising their design). In practice, average performance might be interesting (and computed by an integral). In practice, focusing on worst case behaviour avoids the need to make lots of questionable judgements about relative probability of various conditions. $\endgroup$
    – mikado
    Commented Sep 29, 2022 at 19:46
4
$\begingroup$

The issue with this integrand seems to be the highly oscillatory nature which makes it difficult to find the right step spacing. In the following I will give two methods to deal with the small step sizes.

The integrand:

g[t_] = Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t]

First method: integration using NDSolve

One possibility is to use the adaptive methods of NDSolve to find the right step sizes:

Edit: As recommended in the comments, I changed the starting point from 10^(-12) to $MachineEpsilon/2

AbsoluteTiming[
NDSolveValue[{v'[t] == g[t], v[$MachineEpsilon/2] == 0}, 
v[2 Pi], {t, $MachineEpsilon/2, 2 Pi}]/(6*Pi) - 1]

output: (*{0.391597, -2.92904*10^-6}*)

If we need more precision we can adjust the options of NDSolve, for example :

AbsoluteTiming[
NDSolveValue[{v'[t] == g[t], v[$MachineEpsilon/2] == 0}, 
v[2 Pi], {t, $MachineEpsilon/2, 2 Pi}, PrecisionGoal -> 10, 
AccuracyGoal -> 10, Method -> "ExplicitRungeKutta"]/(6*Pi) - 1] (* relative error *)

output: (*{0.438849, 1.79297*10^-10}*)

Second method: summing over sub-intervals specified by the zeroes of the function.

A more manual approach is to consider a partition of the interval defined by the zeroes of the function. Specifically, for x[r] the r'th zero of the integrand, we consider the following partition of the original interval Interval[0,2*Pi]=IntervalUnion[Interval[x[0],x[1]],Interval[x[1],x[2]],....] (this is just a representation not part of the code)

Find the zeroes:

Using Solve or SolveValues takes a lot of time. As such, Reduce is used instead as it is roughly seven times faster in this case :

a = Reduce[Cos[2022*t]*Sin[10050*t]*Sin[10251*t] == 0 && 0 <= t <= 2*Pi, t];

However, the output format is a mixture of equalities and inequalities for auxiliary constants. The code below extracts and constructs numerical zeroes from the output of Reduce:

zeroes = (List @@ a[[;; -5, 2]])~Join~
  DeleteDuplicates[Apply[Join, (#[[3, 2]] /. C[1] -> Range[0, #[[2, 5]]] &) /@ a[[-4 ;; -1]] ]]

Next, we use Gauss-Legendre points in each sub-interval. The most (computer) efficient way might be to use

NIntegrate`GaussKronrodRuleData[points,precision]

as used in the documentation on integration rules and then rescale the points according to the sub-intervals. However, for brevity and convenience, we will use GaussianQuadratureWeights from the NumericalDifferentialEquationAnalysis` package:

Needs["NumericalDifferentialEquationAnalysis`"]

points = (GaussianQuadratureWeights[21, ##] &) @@@ 
Partition[Sort@zeroes, 2, 1];

pointsf = SortBy[Join @@ points, First];

Total[pointsf[[All, 2]]*g[pointsf[[All, 1]]]]/(6*Pi) - 
1 // AbsoluteTiming (* relative error *)

output: (* {0.151416, 8.65974*10^-13} *)

One should also consider the time used to obtain the zeroes but if the integral has to be used for a range of frequencies, using Reduce symbolically would be a one time cost.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Saw your answer now after posting similar one. You could use v[$MachineEpsilon/2] == 0 instead and this should produce better result. no need to use magic number like 10^(-12) $\endgroup$
    – Nasser
    Commented Sep 28, 2022 at 23:40
  • $\begingroup$ +1. Very interesting and nonstandard. $\endgroup$
    – user64494
    Commented Sep 29, 2022 at 5:33
  • $\begingroup$ @user64494 I might have seen that method applied on stack exchange in the past but I do not remember well. I have used that method in the past and I am still a bit surprised how well it works. Quadrature methods have exponential convergence on analytical functions with a finite radius of convergence. The trapezoidal rule has exponential convergence on periodic functions with a finite radius of convergence and double exponential methods are fairly efficient in handling endpoint singularities. When no reasonably fast exponential convergence is possible, perhaps it is best to use NDSolve. $\endgroup$
    – userrandrand
    Commented Sep 29, 2022 at 12:24
4
$\begingroup$

Stripping "GlobalAdaptive" of the "extras" produces correct results quickly:

AbsoluteTiming[
 res1 =
  NIntegrate[
   Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t], {t, 0, 2*Pi}, PrecisionGoal -> 6, 
   Method -> {"GlobalAdaptive", MinRecursion -> 6, 
     MaxRecursion -> 1000, "MaxErrorIncreases" -> 10^6, 
     "SingularityDepth" -> Infinity, "SymbolicProcessing" -> False}]
 ]

(* {0.630764, 18.8496} *)

6 Pi - res1

(* 2.63789*10^-11 *)

AbsoluteTiming[
 res2 =
  NIntegrate[
   Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t], {t, 0, 2*Pi}, PrecisionGoal -> 12, WorkingPrecision -> 30, 
   Method -> {"GlobalAdaptive", MinRecursion -> 6, 
     MaxRecursion -> 1000, "MaxErrorIncreases" -> 10^6, 
     "SingularityDepth" -> Infinity, "SymbolicProcessing" -> False}]
 ]

(* {7.55792, 18.8495559215387594307749777268} *)

6 Pi - res2

(* 8.825728*10^-22 *)

Let us verify using an alternative method, "LocalAdaptive", and alternative accuracy and precision goals (but still without symbolic pre-processing):

AbsoluteTiming[
 res3 = NIntegrate[
   Cos[2022*t]*Sin[10050*t]*Sin[10251*t]/Sin[50*t]/Sin[51*t], {t, 0, 2*Pi}, AccuracyGoal -> 12, PrecisionGoal -> 12, 
   MaxRecursion -> 10^6, MaxPoints -> 10^8, 
   Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}]]

(* {0.877531, 18.8496} *)

6 Pi - res3

(* -4.57681*10^-9 *)

```
Improve this answer
$\endgroup$
4
  • $\begingroup$ +1 for your extremal options. $\endgroup$
    – user64494
    Commented Sep 29, 2022 at 5:41
  • 1
    $\begingroup$ @AntonAntonov I tested the methods proposed so far and after replacing the upper bound 2*Pi by 10*Pi, this method seems to be the fastest and most precise. $\endgroup$
    – userrandrand
    Commented Sep 29, 2022 at 19:25
  • 2
    $\begingroup$ @userrandrand Thanks. You can post another answer with a table that has the timings. $\endgroup$
    – Anton Antonov
    Commented Sep 29, 2022 at 19:28
  • $\begingroup$ @AntonAntonov sorry it seems that I just had to increase the number of points rather than bisections in the method provided by Goofy. I wrote an answer with the timings but then deleted it as they are not enlightening other than the case of NDSolve which fails starting at 8*Pi for an unknown reason. $\endgroup$
    – userrandrand
    Commented Sep 29, 2022 at 23:06
3
$\begingroup$

It's kind of a pain, but you can break up the integral yourself.

integrand = (Cos[2022*t]*Sin[10050*t]*Sin[10251*t])/(Sin[50*t]*Sin[51*t])

To determine where one might break it up, we can plot pieces of the integrand.

Plot[integrand, {t, 0, .0005}, PlotRange -> All]

enter image description here

and the end has the same behavior.

Plot[integrand, {t, 2 \[Pi] - .0005, 2 \[Pi]}, PlotRange -> All]

enter image description here

Those two plots show a substantial chunk of the integral value.

I broke NIntegrate up as follows:

First the two sections I plotted,

int1 = NIntegrate[integrand, {t, 0, .0005}]
(*5.671497870331018*)

int2 = NIntegrate[integrand, {t, 2 \[Pi] - .0005, 2 \[Pi]}]
(*5.671497870331018*)

The next two segments are 1) the end of the first segment to 0.1 Pi and 2) 1.9 Pi to the start of the last segment.

int3 = NIntegrate[integrand, {t, .0005, .1 \[Pi]}]
(*-0.8981984314029957*)

int4 = NIntegrate[integrand, {t, 1.9 \[Pi], 2 \[Pi] - .0005}]
(*-0.8981984314029957*)

and the final integral goes from 0.1 Pi to 1.9 Pi in 0.1 Pi increments.

int5 = Total@ Table[NIntegrate[integrand, {t, .1 \[Pi] i, .1 \[Pi] i + .1 \[Pi]}], {i, 1, 18}]
(*9.30295704371735*)

Add them all.

int = int1 + int2 + int3 + int4 + int5
(*18.84955592150306*)

Compare to the known value

6 \[Pi] // N
(*18.84955592153876*)

Pretty close. None of the integrations gave me failed to converge errors, so this method evidently works, but I have no idea why NIntegrate cannot do this automatically.

Improve this answer
$\endgroup$
1
  • $\begingroup$ +1. Thank you. This is rather a handwork with help of Mathematica. $\endgroup$
    – user64494
    Commented Sep 28, 2022 at 6:42

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.