# Accuracy of function past 10^(-17)

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?

• For your first NIntegrate based solution did you already try specifying a higher WorkingPrecision? Commented Nov 26, 2018 at 19:12
• You might be interested in this Q&A, which presents several ways, including built-in ones, for computing Gauss-Legendre quadrature. Commented Nov 26, 2018 at 20:00

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.