1
$\begingroup$

I am trying to integrate an analytic function:

Func[Θ_, ϵ_] = (π (28 + (-2 + \
ϵ) ϵ (16 + (-2 + ϵ) ϵ (29 + (-2 + \
ϵ) ϵ (62 + 15 (-2 + ϵ) ϵ))) + 
        2 (-17 + (-2 + ϵ) ϵ (-22 + (-2 + \
ϵ) ϵ (5 + 
                 2 (-2 + ϵ) ϵ (32 + 
                    3 (-2 + ϵ) ϵ)))) Cos[\
Θ] + (-1 + ϵ)^2 (4 + (-2 + ϵ) \
ϵ (12 + (-2 + ϵ) ϵ (59 + 
                 3 (-2 + ϵ) ϵ))) Cos[
          2 Θ] + 
        2 (-1 + ϵ)^4 (1 + 
           4 (-2 + ϵ) ϵ) Cos[3 Θ] - 
        6 (-2 + ϵ)^4 ϵ^4 Csc[Θ/
          2]^2))/(6 (-1 + ϵ)^4 (-1 + (-2 + ϵ) \
ϵ + (-1 + ϵ)^2 Cos[Θ])^2) + 
   Log[(2 - ϵ)/ϵ] (
    4 π (-1 + ϵ)^2 (2 + (-2 + ϵ) ϵ) - 
     2 π (4 + 
        3 (-2 + ϵ) ϵ (2 + (-2 + ϵ) \
ϵ)) Cos[Θ])/(-1 + ϵ)^5 - 
   Log[((1 - ϵ) Sin[Θ/2] + 
      Sqrt[(2 - ϵ) ϵ + (1 - ϵ)^2 Sin[\
Θ/2]^2])/
     Sqrt[(2 - ϵ) ϵ]]/((2 - ϵ) ϵ \
+ (1 - ϵ)^2 Sin[Θ/2]^2)^(5/2) π /(
    4 ((-1 + ϵ)^5) )
     Csc[Θ/
     2]^3 (-(-2 + ϵ)^5 ϵ^5 - (-2 + ϵ)^4 \
ϵ^4 (4 + 5 (-2 + ϵ) ϵ) Sin[Θ/
        2]^2 + 10 (-2 + ϵ)^3 (-1 + ϵ)^2 \
ϵ^3 (-4 + 
         3 (-2 + ϵ) ϵ) Sin[Θ/2]^4 - 
      40 (-2 + ϵ)^2 (-1 + ϵ)^4 ϵ^2 (-6 + \
(-2 + ϵ) ϵ) Sin[Θ/2]^6 + 
      16 (-2 + ϵ) (-1 + ϵ)^6 ϵ (-20 + (-2 \
+ ϵ) ϵ) Sin[Θ/2]^8 + 
      128 (-1 + ϵ)^8 Sin[Θ/2]^10);

I can do it using NIntegrate, but it doesn't give very accurate details:

A = Table[
  1/(2 n + 1)
    NIntegrate[
    Sin[Θ] LegendreP[n, 
      Cos[Θ]] Func[Θ, 1/
      10], {Θ, 0 + 10^-3, π - 10^-3}, 
    PrecisionGoal -> 100, AccuracyGoal -> 100], {n, 1, 50}]

returns

{-1.07674*10^-6, 0.949863, 0.114375, 0.0230668, 0.00586787, \
0.00169616, 0.000530953, 0.000174369, 0.0000597002, 0.0000203903, 
 7.37619*10^-6, 2.21039*10^-6, 8.83215*10^-7, -7.5787*10^-8, 
 3.47186*10^-8, -3.41165*10^-7, -7.25632*10^-8, -3.43366*10^-7, \
-7.9973*10^-8, -3.15529*10^-7, -7.47012*10^-8, -2.8831*10^-7, \
-6.86621*10^-8, -2.6488*10^-7, -6.33192*10^-8, -2.4487*10^-7, \
-5.86535*10^-8, -2.27233*10^-7, -5.43555*10^-8, -2.06886*10^-7, \
-5.18649*10^-8, -1.38998*10^-7, -6.71511*10^-8, -3.63283*10^-7, 
 1.10033*10^-7, -1.21779*10^-6, 1.98552*10^-6, -4.25593*10^-6, 
 5.26844*10^-6, -0.0000319476, 0.000205053, -0.00010609, 0.000150802, \
-0.000543319, 0.000401249, 0.00114871, -0.00146559, 0.0516983, \
-0.00297781, -0.222248}

so somewhere near n = 15, the accuracy can no longer keep up. That's fine, I wrote myself a Gaussian quadrature function to tackle this:

GaussQuadInt[CorrCurve_, ϵ_, 
  Lmin_, Lmax_, nPoints_] := 
 Block[{Coefficients, LegendrePDiff, CorrFunction, GaussLegendreNodes,
    GaussLegendreWeights, CorrFunctionVals},
  LegendrePDiff[l_, x_] = D[LegendreP[l, y], y] /. {y -> x};

  CorrFunction = 
   Compile[{{Θ, _Real}}, 
    Evaluate@CorrCurve[Θ, ϵ]];

  GaussLegendreNodes = 
   Sort[N[x /. Solve[LegendreP[nPoints, x] == 0], 100]];
  GaussLegendreWeights = 
   ParallelTable[
    2/((1 - GaussLegendreNodes[[i]]^2) LegendrePDiff[nPoints, 
      GaussLegendreNodes[[i]]]^2), {i, 1, nPoints}];

  CorrFunctionVals = 
   ParallelTable[
    N[CorrFunction[ArcCos[GaussLegendreNodes[[i]]]], 100], {i, 1, 
     nPoints}];

  Coefficients = Table[
    1/(2 l + 1)
      ParallelSum[
      GaussLegendreWeights[[i]] LegendreP[l, 
        GaussLegendreNodes[[i]]] CorrFunctionVals[[i]], {i, 1, 
       nPoints}]
    , {l, Lmin, Lmax}];

  Return[Coefficients];]

which I then use to perform the same computation:

newCoeffs = GaussQuadInt[Func, N[1/10, 100], 1, 60, 200]

and it gives:

{6.14323*10^-15, 0.949866, 0.114376, 0.0230682, 0.00586816, \
0.00169716, 0.000531168, 0.000175133, 0.0000598703, 0.0000210084, 
 7.51663*10^-6, 2.72966*10^-6, 1.00285*10^-6, 3.71855*10^-7, 
 1.38916*10^-7, 5.22141*10^-8, 1.97254*10^-8, 7.48358*10^-9, 
 2.84942*10^-9, 1.08827*10^-9, 4.16741*10^-10, 1.59952*10^-10, 
 6.15151*10^-11, 2.36992*10^-11, 9.14443*10^-12, 3.53327*10^-12, 
 1.36687*10^-12, 5.29339*10^-13, 2.05214*10^-13, 7.96245*10^-14, 
 3.09582*10^-14, 1.20238*10^-14, 4.69478*10^-15, 1.84874*10^-15, 
 7.27753*10^-16, 2.92705*10^-16, 1.40073*10^-16, 6.39707*10^-17, 
 2.42422*10^-17, 
 4.28862*10^-17, -4.28038*10^-17, -3.89701*10^-17, -5.00079*10^-17, \
-4.22961*10^-17, -3.92505*10^-17, -3.66997*10^-17, -4.47741*10^-17, \
-4.34709*10^-17, -2.59858*10^-17, -4.62965*10^-17, 2.22314*10^-18, 
 3.54379*10^-18, 2.43186*10^-19, 8.99192*10^-19, 
 6.97015*10^-18, -1.23964*10^-18, 1.96099*10^-19, 
 8.13244*10^-18, -3.48402*10^-18, 7.02491*10^-18}

which, unfortunately, loses accuracy too, around n = 40. Here is a plot to illustrate my issue, blue dots are the result of NIntegrate, orange dots are the result of my Gaussian quadrature computation: Fig. 1

What can I do to get a better accuracy of this integration?

$\endgroup$
  • $\begingroup$ For your first NIntegrate based solution did you already try specifying a higher WorkingPrecision? $\endgroup$ – Thies Heidecke Nov 26 '18 at 19:12
  • 1
    $\begingroup$ You might be interested in this Q&A, which presents several ways, including built-in ones, for computing Gauss-Legendre quadrature. $\endgroup$ – Michael E2 Nov 26 '18 at 20:00
3
$\begingroup$

CorrFunction is compiled and thus converts all its inputs into machine precision and also returns numbers in machine precision. All calculations depending on the output of CorrFunction will also be coerced to machine precision.

In order to obtain higher precision results, replace CorrFunction by a conventional function that can take advantage of arbitrary precision numbers. For example with

CorrFunction[θ_] = CorrCurve[θ, ϵ];

I obtain results with Precision close to 95.. And the execution of GaussQuadInt[Func, N[1/10, 100], 1, 60, 200] does not take longer than before.

$\endgroup$
  • $\begingroup$ Great, many thanks. I hadn't realised that compiling the function is what's causing problems. $\endgroup$ – Gordon Nov 26 '18 at 18:37
  • $\begingroup$ You're welcome. $\endgroup$ – Henrik Schumacher Nov 26 '18 at 19:05

Your Answer

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

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