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.341905887336485206837863135597`20.*^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,
Method -> {"GlobalAdaptive",
"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]]}
];
InterpolationOrder -> All
toNDSolve
. This will maketheCentralTrace
more accurate between integration steps; you may not even need such a smallMaxStepSize
. $\endgroup$