### Section 13: Analyzing the Annular Laurent Integrals

$$\newcommand{\bint}{\displaystyle{\int\hspace{-10.4pt}\Large\mathit{8}}} \newcommand{\res}{\displaystyle{\text{Res}}}$$

The reader is asked to review Section 10: Region of convergence of algebraic power series for background on this section.

In this section we study the Laurent integrals of annular Puiseux series: \begin{equation} \begin{aligned} w_n(z)&= \sum_{k=0}^{\infty} a_k (z^{1/n})^k+\sum_{k=1}^{\infty} \frac{b_k}{\left(z^{1/n}\right)^k} \\ &=A(z)+S(z) \end{aligned} \end{equation} with \begin{equation} \begin{aligned} a_k&=\frac{1}{2n\pi i} \bint \frac{w_n(z)}{\left(z^{1/n}\right)^{k+n}} dz\\ b_k&=\frac{1}{2n\pi i} \bint w_n(z)\left(z^{1/n}\right)^{k-n} dz.\\ \label{equation:equation4} \end{aligned} \end{equation} and $A(z)$ being the analytic terms of the series and $S(z)$, the singular terms. Or in symmetrical form: \begin{equation} \begin{aligned} a_k&=\frac{1}{2n\pi i} \bint \frac{w_n(z)}{\left(z^{1/n}\right)^{k+n}} dz\\ w_n(z)&=\sum_{p=-\infty}^\infty a_p (z^{1/n})^p. \end{aligned} \label{equation:equation100} \end{equation}

We use the function in Section 10 \begin{equation} \begin{aligned} f(z,w)&=(-z^2+z^3)\\ &+(-4 z+3 z^2)w\\ &+(-z^3-9 z^4)w^2\\ &+(-2+8 z+4 z^2-4 z^3)w^3\\ &+(6-8 z^2+7 z^3+8 z^4)w^4 \end{aligned} \label{equation:equation20} \end{equation} In Section 10, we computed the branch continuation table for this function. In this section, we analyze the integral expression of (\ref{equation:equation100}) in the following form: \begin{equation} \begin{aligned} c_k&=\frac{1}{2n\pi i} \bint \frac{w_n(t)}{\left(z^{1/n}\right)^{k+n}} dz=\frac{1}{2 n \pi} \int_{t_0}^{t_e} \frac{w_n(t)re^{it}}{\left(re^{it}\right)^{\frac{k+n}{n}}}\\ &=\frac{1}{2 n\pi} \int_{t_0}^{t_e} w_n(t) \left(re^{it}\right)^{-k/n} dt\\ &=\frac{1}{2 n\pi r^{k/n}} \int_{t_0}^{t_e} w_n(t)\big[\cos(tk/n)-i\sin(tk/n)\big] dt\\ &=\frac{1}{2 n\pi r^{k/n}} I(k,n). \end{aligned} \label{eqn5} \end{equation}

We are however immediately confronted with oscillating integrals which increase in frequency as the parameter $k$ is increased. So the question arises: How do we compute (\ref{eqn5}) accurately?

Consider the analytic domains for the branches in the first annulus in the continuation table: \begin{aligned} \text{Branch} & & \text{Domain} \\ \{1,3\} & & (0,0.00915) \\ \{2\} & & (0,0.693) \\ \{4\} & & (0,0.00915) \\ \end{aligned}

And we wish to study the power expansion for the $2$-cycle branch which we can calculate very precisely using Newton Polygon. How will those results compare with the numerical integration of (\ref{eqn5})? Or more importantly, how can we integrate (\ref{eqn5}) in order to achieve a desired accuracy with Newton Polygon as the benchmark so that we can then apply what we learned to the remaining annuli for which we do not have a benchmark?

First consider the integrand for the $\{1,3\}$ branch in annulus $1$. Figure 1 shows the real component of the integrand for various values of $k$ and illustrates how the integrand increases in frequency as $k$ increases.

And to give an example of the issue, consider the default Mathematica calculations of Equation $\ref{eqn5}$ for the $c_1$ coefficient:

Table 1: $c_1$ results for $\{1,3\}$
ItemReal PartImaginary Part
Actual Value01.2053049556196973925392665399622067511862925003555
Computed Value-0.00000009479699058026481.20530489619523
Difference-0.0000000947969905802648-0.00000005942446712481342

and the results for $c_{100}$:

Table 2: $c_{100}$ results for $\{1,3\}$
ItemReal PartImaginary Part
Actual Value0.000000000000000000563533493348688562530
Computed Value-0.0000000009006157096675720.000000002507049845040755
Difference-0.0000000009006157102311050.000000002507049845040755

In the first case we are accurate to $7$ significant figures and in the second case we are off by $9$ digits. Fortunately however, we can significantly improve these results. First though we review the code to implement the integration.

### Setting up the Analysis:

From the continuation table, we know the $\{1,3\}$ branch in annulus $1$ has a domain of $(0,0.0092)$. Our first objective is to select a value of $r$ in which to perform the integration. We can start by simply selecting the mean of the domain or about $0.00451$. However, we cannot simply choose a decimal representation for this radius since Mathematica will treat the value as a machine-precision number with about $15$ digits of precision. This will constrain our integration results to the same level of precision. We can represent the radius as a arbitrary-precision number with $r_0=451/100000$.

Our next objective is to obtain the trace of the branch along a circular route $z(t)=re^{it}$. We do this by solving the monodromy differential equation: $$\frac{dw}{dt}=-\frac{f_z}{f_w}\frac{dz}{dt};\quad w(z_0)=w_0$$ keeping in mind the continuation table was constructed relative to the annular reference points for each annulus. For example, branch $\{1,3\}$ in annulus one begins at $r_0 e^{it}$ with $r_0$ equal to the mean of the domain and $t=\pi/6$ for this example. Therefore, for this branch, $z_0=r_0 e^{\pi/6 i}$. And we can compute $w_0$ by first solving for the roots of $f(z_0,w)=0$, sorting the values first by real part than imaginary part, then selecting from the sorted list, the first value since the branch designation begins with $1$. Having $z_0$ and $w_0$, we then numerically solve the monodromy differential equation over a $4\pi$ route along the path $z(t)=r_0e^{it}$ which in this discussion is called the Central Trace of the branch. However, we must solve for both $w_0$ and the Central Trace at a sufficiently high precision in order to integrate Equation $\ref{eqn5}$ accurately. Remember, Mathematica will perform a calculation at the lowest precision of the input data. In practice we obtain $400$ digits of precision for the root $w_0$. Once we obtain the Central Trace, we then can solve $\ref{eqn5}$.

### Integrating the Annular Laurent Integrals:

Equation $\ref{eqn5}$ is sensitive to several parameters. These include oscillatory behavior, radius of integration, working precision, and integration strategy and rule. In this section we discuss these parameters as we attempt to arrive at an efficient method to accurately evaluate the integrals.

We have two main concerns: (1) Computing an analytically-continuous version of the branch trace, and (2), numerically integrating Equation $\ref{eqn5}$

1. ### Computing the branch trace

We numerically solve the monodromy differential equation using the built-in function $\texttt{NDSolve}$. NDSolve works by taking a sequence of steps and uses an adaptive procedure to determine the size of the steps. In general, if the solution appears to be varying rapidly at a particular spot, NDSolve reduces the step size to better track the solution. The precision or accuracy of the result can be specified and NDSolve will reduce the step size until the solution reached satisfies either the AccuracyGoal or the PrecisionGoal or the number of adjustments exceeds either a default limit or user-specified limit. Alternatively, a MaxStepSize can be specified. NDSolve also uses WorkingPrecision to determine the precision to use in its internal computations. A larger AccuracyGoal or PrecisionGoal generally requires a larger WorkingPrecision. Various methods can be used by NDSolve such as ExplicitEuler or RungeKutta. In this section, we will use the method "StiffnessSwitching" as it was shown through trial-and-error to give the best results.

We will use the following code to analyze the integrals.

#### Mathematica code

(* Mathematica code to compute annular Laurent integral for c_100 for \
branch {1,3} in annulus 1 for the following function*)

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);
(* choose coefficient 100: *)
k = 100;
(* set precision values as needed for precision integration *)
workingPrecision = 30;
maxStepSize = 1/20000;
(* set branch size to 2 as this is a 2-cycle branch *)
theBranchSize = 2;
theExp = k/theBranchSize;
(* compute high-precision value of integral from Newton Polygon *)
c100Val =
3.269654303316744493023319752352130734395351218030940654901174166974\
95400631911260.*^97;
c100Integral = 2 Pi theBranchSize rnorm^theExp c100Val;
(* set radius of annular reference point to a arbibrary-precision \
value *)
rnorm = 0.00459986;
rnorm = Round[rnorm, 10^-10];
(* set up integration limits *)
tStart = Pi/6;
tEnd = tStart + 2 \[Pi] theBranchSize;
(* set up z_0 *)
zstart = rnorm Exp[I tStart];
zend = rnorm Exp[I tEnd];
(* get w_0 for monodromy differential equation *)
theBaseValues =
w /. NSolve[theFunction == 0 /. z -> zstart, w,
WorkingPrecision -> 400];
theBaseValues = Sort[theBaseValues, If[Re[#1] != Re[#2],
Re[#1] < Re[#2]
,
Im[#1] < Im[#2]
] &];

wstart = theBaseValues[];
(* construct monodromy differential equation *)
tDeriv = w'[
t] == ((-(D[theFunction, z]/D[theFunction, w]) (I rnorm Exp[
I t])) /. {w -> w[t], z -> rnorm Exp[I t]});
(* set up block to run analysis with default setting or with higher \
precision *)
precisionFlag = True;
integrationFlag = False;
printFlag = False;
precisionArray = {30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80};
stepSizeArray = {1/5000, 1/5000, 1/5000, 1/5000, 1/5000, 1/5000,
1/5000, 1/5000, 1/5000, 1/5000, 1/20000};

myStatTable3 = Table[
Print["analyzing: ", precisionArray[[index]]];
(* Solve DE *)
difTiming = AbsoluteTiming[
First[
NDSolve[{tDeriv, w[tStart] == wstart}, w, {t, tStart, tEnd},
MaxSteps -> 20000000, InterpolationOrder -> All,
AccuracyGoal -> precisionArray[[index]],
MaxStepSize -> stepSizeArray[[index]],
WorkingPrecision -> precisionArray[[index]],
Method -> "StiffnessSwitching"]]];
theSolution = difTiming[];

theUpperCentralTrace[t_] := Evaluate[Flatten[w[t] /. theSolution]];
theBaseValues =
w /. NSolve[theFunction == 0 /. z -> zend, w,
WorkingPrecision -> 400];
theBaseValues = Sort[theBaseValues, If[Re[#1] != Re[#2],
Re[#1] < Re[#2]
,
Im[#1] < Im[#2]
] &];
diff = theBaseValues[] - theUpperCentralTrace[tEnd];
Print["Dif: ", diff];
(* solve integral *)
theIntegralTiming =
AbsoluteTiming[
NIntegrate[
theUpperCentralTrace[t] (Cos[t theExp] - I Sin[t theExp]), {t,
tStart, tEnd},
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 400,
Method -> "LevinRule"},
WorkingPrecision -> precisionArray[[index]],
AccuracyGoal -> precisionArray[[index]],
MaxRecursion -> Infinity]];
intResults = theIntegralTiming[];
intDif = c100Integral - intResults;
Print["int dif: ", intDif];
(* store results *)
{difTiming[], theIntegralTiming[], precisionArray[[index]],
stepSizeArray[[index]], theUpperCentralTrace[tEnd], diff,
intResults, intDif}, {index, Length@precisionArray}];



We first solve the differential equation using increasing precision and then compare the actual value of $w(\text{tEnd})$ at the end of the calculation with the value calculated by NDSolve. It is known that the error of NDSolve increases as $t$ increases so we use the value at the end of the calculation as an upper estimate of the error. Figure 2 displays the results of these calculations as PrecisionGoal is raised from $30$ to $80$. In the code above, a maximum step size is specified as MaxStepSize=1/5000 except for the last run where MaxStepSize=1/20000. In the diagram for example, when PrecisionGoal is at $40$, it takes about $5$ minutes for NDSolve to run and the resulting calculation is accurate to around $30$ decimal digits. With PrecisionGoal set to $80$ and StepSize=$1/20000$, NDSolve took about $40$ minutes to execute and resulted in a value of $w(\text{tEnd})$ accurate to $55$ decimal digits.

To get a better perspective of the accuracy of the NDSolve calculations, Figure 3 tabulates the numerical value of the real component of each run compared to the more precise value returned by Newton Polygon. The Newton Polygon result is shown on the top line of the tables below. The values are compared digit by digit with red digits representing the error. The precision of each run is in the first column of the figure.

$$\begin{array}{l} & -0.0331428140148972385092410337265328053748037369454951841916394\\ 30 & -0.03314281401489723850924103372\color{red}{7}\\[-1ex] 35 &-0.0331428140148972385092410337265\color{red}{4639}\\[-1ex] 40 &-0.0331428140148972385092410337265\color{red}{404417096}\\[-1ex] 45 &-0.0331428140148972385092410337\color{red}{56895089729476667}\\[-1ex] 50 & -0.03314281401489723850924103372653280\color{red}{454408519388436}\\[-1ex] 55 &-0.033142814014897238509241033726532805374\color{red}{7927933023310298}\\[-1ex] 60 &-0.033142814014897238509241033726532805374\color{red}{797699394890660185003}\\[-1ex] 65 &-0.033142814014897238509241033726532805374\color{red}{8040787735710194028485}\\[-1ex] 70 &-0.03314281401489723850924103372653280537480373694\color{red}{98513142233811}\\[-1ex] 75 &-0.03314281401489723850924103372653280537480373694549\color{red}{64512623050}\\[-1ex] 80 &-0.03314281401489723850924103372653280537480373694549518\color{red}{72940900}\\[-1ex] \end{array}$$
Figure 3: Accuracy of NDSolve at $w(25\pi/6)$

2. ### Computing the annular Laurent integral

In this section we use the built-in numerical integrator $\texttt{NIntegrate}$ to evaluate $\ref{eqn5}$. The primary features of NIntegrate are integration strategy and rule. An integration strategy is an algorithm that prescribes how to improve the integral estimate by creating additional sample points. In the code above, we use GlobalAdaptive. A global adaptive strategy reaches the required precision and accuracy goals of the integral estimate by recursive bisection of the subregion with the largest error estimate into two halves, and computes integral and error estimates for each half. This bisection is continued until the precision or accuracy goal is reached or a maximum number of bisections is reached.

An integration rule computes an estimate of an integral over a region using a weighted sum and provides both an integral estimate and an error estimate as a measure of the integral estimate's accuracy. For highly-oscillatory integrands, the Levin rule is often used. A Levin rule estimates the integral of an oscillatory function by approximating the antiderivative as a product of a polynomial and the oscillatory part of the integrand. This rule is used in the code above.

After we obtain the central trace of the branch at a particular precision, we then solve the integral $$\int_{\pi/6}^{25 \pi/6} \text{theCentralTrace[t]}\big[\cos(100/2 t)-i\sin(100/2 t)\big] dt$$ using the following code:

#### Mathematica code

NIntegrate[theUpperCentralTrace[t] (Cos[t theExp] - I Sin[t theExp]), {t,tStart, tEnd},
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 400,Method -> "LevinRule"},
WorkingPrecision -> precisionArray[[index]],
AccuracyGoal -> precisionArray[[index]],
MaxRecursion -> Infinity]];
`
with precisionArray set to a set of precision values we wish to solve the integral. We then analyze the integrand and combine the performance data for both NIntegrate and NDSolve shown in Figure 3. In the figure, the timing for NIntegrate is less than one minute so the purple line is shown along the precision axis. The figure shows the integration precision closely follows the precision of the NDSolve results but takes much less time: NDSolve at a precision of 80 took 40 minutes to compute but the integral at the same precision took about 11 seconds.

Returning now to the results in Table 2 above, we re-calculate the integral for $c_{100}$ with PrecisionGoal set to the values in Figures $4$ and $5$. These show a significant improvement in the results albeit at the expense of greater execution time; the NDSolve calculation at PrecisionGoal set to 100 with step size of 1/20000 took 70 minutes on a 2.2 GHz machine. A goal of 110 with step size of 1/40000 took almost four hours.

$$\begin{array}{l} & 0.000000000000000000563562896773434577333144822440739446643144430190316760651438 \\ 80 & 0.00000000000000000056356289677343457733314482244073944664\color{red}{2945459458066021398616}\\[-1ex] 100 & 0.000000000000000000563562896773434577333144822440739446643144430190\color{red}{238024423890}\\[-1ex] 110 & 0.000000000000000000563562896773434577333144822440739446643144430190316760\color{red}{495288} \end{array}$$
Figure 4: Accuracy of Laurent intergral for $c_{100}$

Now we compute the value of $c_{100}$ from Equation $\ref{eqn5}$ using the integral results and compare it to the results computed by Newton Polygon:

$$\begin{array}{l} & 32696543033167444930233197523521307343953512180309406549011741669749540063191123210966809182168046.0\\ 80 & 326965430331674449302331975235213073439\color{Red}{41968384316760570075700000000000000000000000000000000000000.}\\[-1ex] 100 & 3269654303316744493023319752352130734395351218030\color{Red}{4838465365700000000000000000000000000000000000000.}\\[-1ex] 110 & 326965430331674449302331975235213073439535121803094065\color{red}{39952300000000000000000000000000000000000000.} \end{array}$$
Figure 5: Accuracy of $c_{100}$

Even at a precision of 110 digits, the integral value is still quite a bit off. However, this is a worst-case example and our primary goal here is not so much to compute the power expansions extremely accurately, but rather study the accuracy of the numerical integrations and how the accuracy is affected by step size and working precision. We will find later that even though the integral results above for $c_{100}$ was much different than the expected value, we can still obtain very precise results for both the Root Test and power expansion. For example, using the values in Figures 4 and 5, we compute the value of the $1/50$'th root for both expressions using $1/\left(|c_{100}|\right)^{1/50}$ for the Newton Polygon result and the expression $\text{rnorm}\left(\frac{4\pi}{I(110,2)}\right)^{1/50}$ for the Laurent integral result (see link above for a derivation of these expressions) to obtain

$$\begin{array}{l} & 0.01121269382044745210777383053090085626286627114459173120072541256168416953266 \\[-1ex] 110 & 0.0112126938204474521077738305309008562628662711445917312007\color{red}{875480429857440547659588663698} \end{array}$$
which is much more precise than we need since we will use less than 5 or 6 digits past the decimal place for our work with the Root Test. And in the case of the actual value of $c_{100} z^{50}$ with $z=4/1000 e^{\pi i/4}$ for example, the difference between the Newton Polygon results and the integral results is on the order of $10^{-10}$ with the difference between the actual branch value at this point and the power series computed above being less than $10^{-8}$. These and other results using \ref{eqn5} will be further explored in Section 14.