# Can I obtain more precise integration of highly oscillating integral?

I would like to numerically integrate precisely, $$\frac{1}{2 \pi r^{n}}\int_0^{2\pi} e^{-n i t}C(t)dt$$ with $$C(t)$$ a single-cycled branch of the function $$w(z)$$ given by the implicit expression $$-z^2+z^3+w (-4 z+3 z^2)+w^3 (-2+8 z+4 z^2-4 z^3)+w^2 (-z^3-9 z^4)+w^4 (6-8 z^2+7 z^3+8 z^4)=0$$ and $$n=100$$. $$C(t)$$ is an NDSolve result in the code below.

A plot of Re[e^{100 it] C[t]] is shown below.

I know the precise value of the integral by another method. The following code uses LevinRule and can only get the results accurate to $$10^3$$ difference when I set the working precision to 70 and MaxStepSize to 1/60000. If I attempt to increase working precision or decrease step size, the results are less accurate. I was wondering if there is a way to get the results accurate to less than 1?

Here is an example output for getIntegral[20,20,1/1000] routine listed below for n=20 with working precision 20 and max step size=1/1000: $$\begin{array}{ccccc} 20 & \{20,1/1000\} & \begin{array}{ll} (1.4347) \\ (1.43468)\end{array} & 1.20279* 10^{-8} & 1.52 \end{array}$$ Reading from left to right, n=20, then the working precision and step size, then next column is the integral results below which is the actual value, then next column is the difference and last column is time in seconds.

For

getIntegral[100,70,1/60000],


the results are: $$\begin{array}{ccccc} 100 & \{70,1/60000\} & \begin{array}{ll} (834190592984+17i) \\(834190588733) \end{array} & 4.2 * 10^{3} & 250 \end{array}$$

Note: the result is very large because r=1/2 and recall AccountingForm uses parenthesis to note negative numbers.

theFunction = -z^2 + z^3 + w (-4 z + 3 z^2) +
w^3 (-2 + 8 z + 4 z^2 - 4 z^3) + w^2 (-z^3
- 9 z^4) +w^4 (6 - 8 z^2 + 7 z^3 + 8 z^4)
theBaseValues =
w /. NSolve[theFunction == 0 /. z ->
zstart, w,WorkingPrecision -> 200];
theBaseValues = Sort[theBaseValues,
If[Re[#1] != Re[#2],
Re[#1] < Re[#2]
,
Im[#1] < Im[#2]
] &];
wstart = theBaseValues[[3]];
rnorm = 1/2;
tStart = 0;
tEnd = 2 \[Pi];
zstart = rnorm Exp[I tStart];
actualValue =
-8.34190588733648520683786313559720.*^11;

wDeriv = w'[t] == ((-(D[theFunction, z]/
D[theFunction, w]) (I rnorm Exp[I t]))
/. {w -> w[t],z -> rnorm Exp[I t]});

getIntegral[j_, wp_, sSize_] :=
Module[{numIndex, difResults, myazsol,
theCentralTrace, n1,finalValue,
intResults},

myazsol =
First[NDSolve[{wDeriv, w[tStart] ==
wstart}, w, {t, tStart, tEnd},
MaxSteps -> 20000000, MaxStepSize ->
sSize,WorkingPrecision -> wp]];
theCentralTrace[t_] = Evaluate[Flatten[w[t]
/. myazsol]];
n1 = AbsoluteTiming[
NIntegrate[( Exp[I t])^-j
theCentralTrace[t], {t, tStart, tEnd},
WorkingPrecision -> wp,
"MaxErrorIncreases" -> 10000,
Method -> "LevinRule"}, MaxRecursion ->
500]];
finalValue = 1/(2  \[Pi] rnorm^j) (n1[[2]]);
intResults = actualValue;
difResults = Abs[intResults - finalValue];
{j, {wp, sSize},
Column[{N[AccountingForm[finalValue], 20],
N[AccountingForm[intResults], 20]},
Alignment -> Left],
N[ScientificForm[difResults], 10], n1[[1]]}
];


• I suspect this problem will have the same solution that my earlier question on a related system did: mathematica.stackexchange.com/questions/186829/… – Kevin Ausman Apr 14 at 14:55
• I suspect I could achieve greater precision if I were able to raise the WorkingPrecision of NDSolve. However, when I attempt to do this to say 105, the integration fails. – Dominic Apr 14 at 17:33
• @Dominic Try adding InterpolationOrder -> All to NDSolve. This will make theCentralTrace more accurate between integration steps; you may not even need such a small MaxStepSize. – Michael E2 Apr 14 at 18:21

The following gives a result that agrees with all the digits for the j = 100 example.

NDSolve seems to be OK if you add the InterpolationOrder -> All option. The "InterpolationPointsSubdivision" preprocessor of NIntegrate seems like the right move, but it was much faster and with a much better result to do a simple, straightforward implementation of an integration rule (losing the protection of any of the checking the NIntegrate would perform).

getIntegral[j_, wp_, sSize_] :=
Module[{numIndex, difResults, myazsol, theCentralTrace, n1,
finalValue, intResults, nodes, weights, errweights},
PrintTemporary["NDSolve"];    (* for impatient people like me :) *)
myazsol =
First[NDSolve[{wDeriv, w[tStart] == wstart}, w, {t, tStart, tEnd},
MaxSteps -> 20000000, MaxStepSize -> sSize,
WorkingPrecision -> wp, InterpolationOrder -> All]];         (* N.B. *)
theCentralTrace[t_] = Evaluate[Flatten[w[t] /. myazsol]];
PrintTemporary["Integrate"];  (* for impatient people like me :) *)
(* Pick a favorite integration rule of sufficiently high order *)
{nodes, weights, errweights} = NIntegrateGaussRuleData[7, wp];
n1 = AbsoluteTiming[
Total[  (* totals the integrals over subintervals *)
Block[{t = Rescale[nodes, {0, 1}, #],                         (* scale nodes to subinterval # = (x_{i-1}, x{i}) *)
$$MinPrecision = wp,$$MaxExtraPrecision = wp},             (* keeps working precision at wp during the dot product *)
weights.(Exp[I t]^-j theCentralTrace[t])*(#[[2]] - #[[1]]) (* sum w_j f(t_j) * (x_{i} - x_{i-1}) = integral over (x_{i-1}, x{i}) *)
] & /@ Partition[Flatten@theCentralTrace["Grid"], 2, 1]    (* maps (/@) integration code over the intervals created by Partition[] *)
(* theCentralTrace["Grid"] yields the steps {{x_0}, {x_1},...} *)
]];                                                           (* Partition transforms {x_0, x_1, x_2,...} to {{x_0, x_1}, {x_1, x_2},...} *)
finalValue = 1/(2 \[Pi] rnorm^j) (n1[[2]]);
intResults = actualValue;
difResults = Abs[intResults - finalValue];
{j,
{wp, sSize},
Column[{
N[AccountingForm[finalValue], 20],
N[AccountingForm[intResults], 20]}, Alignment -> Left
],
N[ScientificForm[difResults], 10],
n1[[1]]}];

PrintTemporary@Dynamic@{Clock[Infinity]}; (* for impatient people like me :) *)
getIntegral[100, 70, 1/5000]


If you want to measure the error estimate from an appropriate integration rule, makes changes like these:

{nodes, weights, errweights} = NIntegrateGaussKronrodRuleData[5, 70];

n1 = AbsoluteTiming[
Total[
Block[{t = Rescale[nodes, {0, 1}, #],
$$MinPrecision = 70,$$MaxExtraPrecision = wp},
{weights, errweights}.(Exp[I t]^-j glurg[t])*(#[[2]] - #[[1]])
] & /@ Partition[Flatten@theCentralTrace["Grid"], 2, 1],
Method -> "CompensatedSummation"
]];
finalValue = 1/(2 \[Pi] rnorm^j) (n1[[2, 1]]);


Update: Combining @Roman's W version of theCentralTrace, we get a speedy solution (the 3*j nodes is heuristic -- might want somewhat higher):

{nodes, weights, errweights} = NIntegrateGaussRuleData[3*j, wp];
n1 = AbsoluteTiming[
Block[{t = Rescale[nodes, {0, 1}, {0, 2 Pi}],
$$MinPrecision = wp,$$MaxExtraPrecision = wp},
weights.(Exp[I t]^-j (W /@ (rnorm Exp[I t])))*(2 Pi)
]];


And then:

getIntegral[100, 60, 1/3000]


• I assume the actualValue for j = 20 has not been given. The above reproduces most of that example, but not the actual value or a believable error estimate. – Michael E2 Apr 14 at 20:19
• Thanks guys! You would not believe how long I have worked on this myself with barely acceptable results. Really did not think it possible to obtain such precise results as Michael's code and do so quickly. Will really help me move forward in my work with algebraic functions and their power expansions. – Dominic Apr 14 at 21:36
• @Dominic I added some notes on the integration code, which might be helpful. – Michael E2 Apr 15 at 15:55
• Thanks Michael. The comments help a lot. – Dominic Apr 15 at 18:11

A much easier way of solving the equation for $$w(z)$$ on the correct branch, instead of integrating the differential equation, is

W[z_?NumericQ] := First@MinimalBy[w /. NSolve[
(-1+z)z^2+w*z(-4+3z)-w^2*z^3(1+9z)+w^3(-2+4z(2+z-z^2))+w^4(6+z^2(-8+z(7+8z)))==0,
w], Abs, 1]


Then you can get the $$n=20$$ solution simply with

With[{n = 20, r = 1/2},
1/(2π r^n) NIntegrate[E^(-I n t) W[r E^(I t)], {t, 0, 2π}]]


-1.43468 - 4.27557*10^-11 I

The $$n=100$$ integral is still difficult, and @MichaelE2 's solution can probably help.

Are you insisting on doing this integral numerically, or are you just interested in the results? If it's the latter, then you can easily get exact results for all values of $$n$$ with a series-expansion of the Root object R[z] below that formally describes the solution of the polynomial equation. It's just a matter of picking the right branch of the root (here, the first branch):

eq[w_, z_] = (-1+z)z^2+w*z(-4+3z)-w^2*z^3(1+9z)+w^3(-2+4z(2+z-z^2))+w^4(6+z^2(-8+z(7+8z)));
R[z_] = Root[eq[#, z] &, 1];
CoefficientList[Series[R[z], {z, 0, 100}], z]


{0, -1/4, 9/128, 85/4096, 131/32768, ...}

Just for $$n=20$$:

SeriesCoefficient[R[z], {z, 0, 20}]


$$-\frac{56833546863764806539901668529}{39614081257132168796771975168}$$

% // N


-1.43468

Just for $$n=100$$:

SeriesCoefficient[R[z], {z, 0, 100}]


$$-\frac{21333059674656860988913490423529488435269036857164193142278059654155\ 9662040070438385095120975342868607822726196631896648747423125210337484\ 60127156003492355000697}{2557336412418860835947804450646561837669251598471144366783821381325\ 1045284411519960025547596296126227741302219746563054759509816764729633\ 229129121792}$$

% // N


-8.34191*10^11

Even for $$n=1000$$ this works with a bit of patience: (not showing the exact result because it's huge)

SeriesCoefficient[R[z], {z, 0, 1000}] // N


-6.7183*10^153

# Update: Version 12

With the new Version 12 function AsymptoticSolve we can get the series-expansion directly, without going through a Root object and without branch warnings. Assuming that you know you want the real-valued branch that goes through zero,

eq[w_, z_] = (-1+z)z^2+w*z(-4+3z)-w^2*z^3(1+9z)+w^3(-2+4z(2+z-z^2))+w^4(6+z^2(-8+z(7+8z)));
AsymptoticSolve[eq[w, z] == 0, {w, 0}, {z, 0, 10}, Reals]


{{w -> -z/4 + 9 z^2/128 + 85 z^3/4096 + 131 z^4/32768 - 444991 z^5/4194304 + 1642905 z^6/134217728 - 19786821 z^7/1073741824 - 662463103 z^8/17179869184 - 542911401095 z^9/4398046511104 - 4656077458645 z^10/140737488355328}}

The list of series coefficients is

CoefficientList[w /.
First@AsymptoticSolve[eq[w, z] == 0, {w, 0}, {z, 0, 10}, Reals], z]


{0, -1/4, 9/128, 85/4096, 131/32768, -444991/4194304, 1642905/134217728, -19786821/1073741824, -662463103/17179869184, -542911401095/4398046511104, -4656077458645/140737488355328}

This is about ten times faster than the previous method using Root & Series.

• I put ours together in my latest edit. (+1) – Michael E2 Apr 14 at 22:21
• That is outstanding Roman! However a lot of my work deals with series expansions of annular branches of the function so I am left with computing these singular series via numerical integration. So I first wanted an efficient algorithm I could check with known results, expansions about the origin, then apply that numerical method to the annular series for which I do not know in advance what the answer is. – Dominic Apr 15 at 11:01
• @Dominic Maybe if you reformulate your question in terms of "How can I do an annular series expansion of this particular function?" someone may be able to help you get around these difficult numerical integrals. Particular the new version-12 function AsymptoticSolve looks promising: reference.wolfram.com/language/ref/AsymptoticSolve.html – Roman Apr 15 at 11:07
• @Roman, thanks for that. Will look into it. However yours and Michael's code is quite complex for me so will first study every detail of it. Will take a while I suspect. You guys are great! – Dominic Apr 15 at 11:36
• The command: CoefficientList[Series[R[z], {z, 0, 100}], z] using all four roots only picks out one single cycle branch and the 2-cycle branch; the 1 and 2 roots generate the same single-cycle series. It should generate two series for the two single-cycle branches and then two more (different ones) for the 2-cycle branch. Can check with Table[Series[Root[eq[#, z] &, n], {z, 0, 5}], {n, 1, 4}]. Looks like a bug. – Dominic Apr 15 at 13:01