Accuracy of function past 10^(-17)

Question

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:

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

Answer 1

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.