Work in progress . . .
This page reviews algorithms and Mathematica code used to determine the radii of convergence of power expansions around singular points of an algebraic function $w(z)$ defined implicitly by the irreducible polynomial expression: $$ \begin{equation} f(z,w)=a_0(z)+a_1(z)w+a_2(z)w^2+\cdots+a_n(z)w^n=0 \end{equation} $$ with $z$ and $w$ complex variables and $a_i(z)$ polynomials with rational coefficients. Details are described in paper Determining the radii of convergence of power expansions around singular points of algebraic functions Readers are encouraged to review the paper first.
Brief overview
Power expansions of $w(z)$ are usually computed by Newton Polygon. Since Mathematica ver 12.1, the builtin function $\texttt{AsymptoticSolve}$ can compute the expansions although there are some issues. In order to determine the radius of convergence of a power expansion, convergencelimiting singular points (CLSPs) are identified for each branch sheet over a base singular point $s_b$ by the method of analytic continuation: the sheets are numericallyintegrated across successively distant singular points until each branch sheet analytically continues onto either a multivalued branch sheet or a single valued sheet with a pole. The singular point where this occurs is the CLSP for the branch sheet. However care must be taken since branch sheets of multicycle branches can have different CLSPs. The first CLSP which this occurs is the CLSP for the branch and establishes the radius of convergence of its associated power expansions.
Steps used to determine and check CLSPs:
 Compute the set $\{s_i\}$ of singular points of $w(z)$. For each $s_i$, sort remaining singular points in order of increasing distance from $s_i$.
 Compute Puiseux series centered at each singular point to a desired working precision and identify branch cycle types according to the categories described in the reference.
 Construct analyticallycontinuous routes between a base singular point,$s_b$, domain and remaining singular point domains.
 Numerically integrate $w(z)$ from $s_b$ domain to $s_i$ domain in order of increasing distance from $s_b$.
 Identify branch sheets over the base singular point, $s_b$ that are analyticallycontinuous across successive singular points over analytic $1$cycle branches until all branch sheets have encountered convergencelimiting singular points (CLSP).
 Apply Root Test to check CLSPs.
We will use the following example from the paper: $$ \begin{equation} f_1(z,w)=z^3 + (z + z^2 + z^3) w + 2 z^2 w^2 + (1 + z + z^3) w^3=0 \end{equation} $$

Computing singular points,singular regions and singular reports
The following code uses $\texttt{NSolve}$ to compute the singular points of $f_1$ then sorts the singular points in order of increasing distance from the origin, computes the $1/3r$ singular regions and then generates a report of the singular reports with poles highlighted in red and a plot in the $z$plane showing each singular point and its $1/3r$ singular region with $r$ the distance to the nearest singular point. Note the singular points are computed to $800$ digits of precision, and since the singular points at $0.68$ and $0.75$ are very close to one another, the associated $1/3r$ regions are very small so the singular regions for these points in the plot below are overwritten by the black dots.
Mathematica code 1
defaultPrecision=800; theFunction=z^3+(z+z^2+z^3) w+2z^2 w^2+(1+z+z^3) w^3; (* compute singular points and sort relative to origin *) theSingularSet=Sort[z/.NSolve[Resultant[theFunction,D[theFunction,w],w]==0, z,WorkingPrecision>defaultPrecision],Abs@#1<= Abs@#2&]; theSingularList=DeleteDuplicates[theSingularSet,Abs[#1#2]<10^50&]; (* find poles *) wExpList=Exponent[theFunction,w,List] poleFunction=Coefficient[theFunction,w,Max@wExpList] thePoleList=z/.NSolve[poleFunction==0,z,WorkingPrecision>defaultPrecision]; polePositions=Flatten[(Position[theSingularList,x_/;Abs[x#]<10^8]&/@ thePoleList)]; (* generate singular list report with red poles *) singularReport=theSingularList; theArgList=Map[Arg,theSingularList]; For[i=1,i<=Length@singularReport,i++, If[MemberQ[polePositions,i], singularReport[[i]]=Style[singularReport[[i]],Red,Bold]; ]; ]; (* compute 1/3r region around each singular point *) nearestSing=(NearestTo[#1,2][theSingularList]&/@ theSingularList)[[All,2]]; nearestDist=MapThread[1/3Abs[#1#2]&,{theSingularList,nearestSing}]; For[i=1,i<=Length@nearestDist,i++, If[nearestDist[[i]]>10, nearestDist[[i]]=10; ]; ]; minRadiusArray=nearestDist; (* generate graphics for singular regions *) minRadiiG=Graphics@MapThread[Circle[ReIm@#1,#2]&,{theSingularList,minRadiusArray}]; (* generate singularListReport *) singularListReport=Row[{Range[Length@theSingularList]//MatrixForm,N[singularReport,6]//MatrixForm}] theSingularListG=Graphics@{PointSize[0.0105],Black,Point@(ReIm@#&/@ theSingularList)}; (* show graphics of singular points and singular regions *) Show[{minRadiiG,theSingularListG},Axes>True,PlotRange>2]

Computing and sorting Puiseux series centered at singular points
In most cases, we can use $\texttt{AsymptoticSolve}$ to compute Puiseux expansions centered at singular points for moderately complex expressions of $(1)$. An expansion about the origin up to degree $512$ took about $8$ seconds on a $4$ GHz desktop:
Mathematica code 2
(* compute Puiseux series and sort at singular reference point *) sCenter = theSingularList[[1]]; AbsoluteTiming[ pSeries = w /. AsymptoticSolve[theFunction == 0, w, {z, sCenter, 512}]; ] minR = minRadiusArray[[sIndex]]; minPrec = Min[Precision@pSeries, Precision@currentSing]; pSeries = SetPrecision[pSeries, minPrec]; currentSing = SetPrecision[currentSing, minPrec]; pSeries = pSeries /. (z  currentSing) > z; theBaseSeriesFun[z_] = pSeries; theBaseValueList = MapThread[{#1, #2} &, {Range[Length@pSeries], pSeries /. z > minR}]; theSortOrder = Sort[theBaseValueList, If[Re[#1[[2]]] != Re[#2[[2]]], Re[#1[[2]]] < Re[#2[[2]]] , Im[#1[[2]]] < Im[#2[[2]]] ] &]; pSeries = pSeries[[theSortOrder[[All, 1]]]]; Precision@pSeries pSeries[[All, 1 ;; 8]] // N // MatrixForm
However in the case of a finiteprecision singular point as the remaining singular points are, the precision of the series is limited by the precision of the singular points: The singular points above were computed to $800$ digits of precision. For example, computing a $512$ degree expansion around $s_5=0.1961+0.5259i$ with $800$ digits of precision produces series with $369$ digits of precision with the precision of the series decreasing with increasing degree of the expansion.
And its important to realize the ramification of the branching may not appear until a sufficient number of terms are generated. Take for example the series: $$ z+z^2+z^3+z^4+z^5+z^6+z^7+z^8+z^9+z^{10}+z^{21/2}+\cdots $$ Had we generated only $10$ terms, we would have concluded this is a Taylor series for a $1$cycle branch.

Constructing an analyticallycontinuous route from $s_b$ to next nearest singular point $s_i$
We first consider continuations at the origin. And in Table $1$, the singular points were arranged in order of absolute value so are also in order of increasing distance to the singular point at the origin. The nearest singular point(s) to the origin is the set $\{0.38520.2530i,0.3852+0.2530i\}$ and our objective is to create an analyticalllycontinuous contour from the singular region around the origin to the singular region around the first singular point in this set. We first consider the general case shown in Figure $2$.
So that the next step is to devise a function to compute points $A$ and $D$ given two singular points and the $1/3r$ singular regions surrounding them. It's not difficult to construct the equation of a line passing through both singular points in terms of the real and imaginary values via an equation of a line $y_b(x)=m(xx_1)+y_1$. And since we know both the singular values $s_b=x_b+iy_b$ and $s_n=x_n+iy_n$ and radii $r_b$ and $r_n$ of each singular radius, it's a simple matter to construct equations for the perimeter of the singular regions: $$ \begin{align*} (xx_b)^2+(yy_b)^2&=r_b^2 \\ (xx_n)^2+(yy_n)^2&=r_n^2 \end{align*} $$ and then find the intersection of $y_b(x)$ with the two circles and from those values, points $A$ and $D$ in the diagram. Once we obtain $A$ and $D$ we can parameterize a path from $A$ to $D$ via $$ \begin{equation} z(t)=A(1t)+D t,\quad 0\leq t \leq 1 \end{equation} $$ Note: The point $F$ in Figure $2$ is used in the case of integrating over a removable singular point. See cited paper for more information about this case.The Mathematica code page Mathematica Code includes the routine $\texttt{getIntegrationData}$. This function takes a base singular point and next singular point and radii of the regions and returns the points $A$, $D$ and $F$ in Figure $2$ with a desired precision. In the case of $s_b=0$ and $s_2=0.35820.2530i$, both region radii are $r_b=r_2=0.14618$ and were computed above to $800$ digits of precision. In order to compute the integration points $A$,$D$, and $F$, we call: $$ \texttt{{pointA,pointD,pointF}=getIntegrationData[s_b,r_b,s_n,r_n,100]} $$ This returns: $$ \texttt{pointA=0.119405  0.0843378 I}\\ \texttt{pointD=0.238809  0.168676 I}\\ \texttt{pointF=0.477618  0.337351 I}\\ $$ We then can construct the integration path from $A$ to $D$ as: $$ \texttt{z(t)=pointA(1t)+pointD t} $$ providing an integration path accurate to $100$ digits of precision.

Numerically integrating the monodromy differential equation
Once we have the integration path from $s_b$ to the next sequential singular point $s_n$, we next analytically continued each branch sheet from the $s_b$ domain to the $s_n$ domain via numerical integration of the monodromy differential equation: Given $f(z,w)=a_0(z)+a_1(z)w+\cdots+a_n(z)w^n=0$, then $\displaystyle \frac{dw}{dt}=\frac{f_z}{f_w}\frac{dz}{dt}$. A set of initial value problems is next set up, one for each branch sheet with initial values at $A$ for each branch sheet and then integrated along the path shown in Figure $2$ to the point $D$. Now recall the cycle sequence over $s_b$ was $\{2,1,2\}$, that is $P_1$ is a $2$cycle series, $P_2$ is a $1$cycle series and $P_3$ is the second sheet of the $2$ cycle series over $s_b$. We next integrate the following $3$ initial value problems: $$ \begin{equation} \frac{dw}{dt}=\frac{f_z}{f_w}\frac{dz}{dt}, \quad\; w_i(0)=P_i(As_b); \quad i=1,2,3 \end{equation} $$ with $z(t)=\text{pointA}(1t)+\text{pointD} t;\quad 0\leq t\leq 1$ and each $P_i(z)$ is a Puiseux series at $s_b$.
We first compute the Puiseux series over $s_2=0.35820.2530i$ up to degree $512$ and sort:
Mathematica code 3
sCenter = theSingularList[[2]]; AbsoluteTiming[ qSeries = w /. AsymptoticSolve[theFunction == 0, w, {z, sCenter, 512}]; ] Precision@qSeries minR = minRadiusArray[[2]]; minPrec = Min[Precision@qSeries, Precision@sCenter]; qSeries = SetPrecision[qSeries, minPrec]; qSeries = qSeries /. (z  sCenter) > z; theNextSeriesFun[z_] =qSeries; theNextValueList = MapThread[{#1, #2} &, {Range[Length@qSeries], qSeries /. z > minR}]; theSortOrder = Sort[theNextValueList, If[Re[#1[[2]]] != Re[#2[[2]]], Re[#1[[2]]] < Re[#2[[2]]] , Im[#1[[2]]] < Im[#2[[2]]] ] &]; qSeries = qSeries[[theSortOrder[[All, 1]]]]; qSeries[[All, 1 ;; 8]] // N // MatrixForm
Mathematica code 4
zTrace[t_]:=pointA(1t)+pointD t; wDeriv=w'[t]==(((D[theFunction,z]/D[theFunction,w]) D[zTrace[t],t])/.{w>w[t],z>zTrace[t]}); baseTest=baseSeries/.z>baseNDZ; theW=NDSolveValue[{wDeriv,w[0]==#},w,{t,0,1},MaxSteps>2000000,MaxStepSize>1/1000,WorkingPrecision>50]&/@ baseTest; #[1]&/@ theW//N nextSeries/.z>(nextNDZtheSingularList[[2]])//N
This was of course a simple example. If one or more branches had continued over a singular point, we would then check the next nearest singular point in the same way and determine if any of those branches continued across it and if so continue checking singular points until all branches had reached convergence limiting singular points.

Checking results with the Root Test
We can apply the Root Test by including the cycle size: $$ \begin{equation} R=\frac{1}{\displaystyle \liminf_{k\to \infty}\; a_k^{\frac{c}{m_k}}} \label{eqn:eqnss2} \end{equation} $$ where $c$ is the cycle size of the series, and the set $\{m_i\}$ is the set of exponent numerators under a least common denominator. Then the $\liminf$ of each Puiseux series can be approximated by forming the set: $$ S=\bigg\{\left(\frac{1}{m_k},\frac{1}{a_k^{\frac{c}{m_k}}}\right)\bigg\} $$ and numerically extrapolating $\displaystyle \lim_{k\to\infty} S$ using a sufficient number of trailing points of $S$. We then approximate $R$ by fitting a suitable curve to the trailing terms of $S$. In some cases, a large number of terms are needed before the sequences settles into a trend (See paper above for an example of this). In our example above, we determied $R=\lefts_2\right\approx 0.43856$ for each series.
The result of the Root test improves with increasing number of terms. Since this is a simple example of degree three, Mathematica can compute a large number of terms in a reasonable amount of time. We therefore will compute the series at the origin to degree $4000$. This means there will be $4000$ terms for the $1$cycle and $8000$ terms for the $2$cycle.
We first compute Puiseux series at the origin up to degree $4000$ and then construct $S$ sets using the following code:
Mathematica code 5
theBranchNum = 3; cycleType = 2; endTerms = 500; branchSeriesF[z_] = pSeries[[theBranchNum]]; seriesLength = Length@branchSeriesF[u]; theexp = Exponent[branchSeriesF[u], u, List]; lcm = LCM @@ (Denominator@theexp); numVals = lcm theexp; clist = Coefficient[branchSeriesF[z], z, theexp]; cTable = Table[ {1/numVals[[k]], N[1/(Abs[clist[[k]]]^(cycleType/numVals[[k]])), 25]}, {k, 5, Length@clist}]; lp1 = ListPlot[cTable, PlotStyle > {PointSize[0.005]}];
This code produces an $S$ set plot shown as the blue scatter points in Figure $3$. Once we have $S$ for each series, its a simple matter to select the lower set of points and fit a curve to the trailing set of $S$. In the code below, routine $\texttt{getLimLine[ctable]}$ select these points, and we then fit a quadratic curve to points $500$ to the last point:
Mathematica code 6
(* code to select lower set of S *) getLimLine[cTable_] := Module[{descentFlag, newList, truncTable, i, oldVal, ascentFlag}, oldVal = cTable[[1, 2]]; newList = {cTable[[1]]}; For[i = 2, i <= Length@cTable, i++, If[cTable[[i, 2]] < oldVal, AppendTo[newList, cTable[[i]]]; oldVal = cTable[[i, 2]]; ]; ]; newList ]; (* now fit data, generate plot and % error report *) newList = getLimLine[cTable[[endTerms ;;]]]; If[Length@newList < 25, Print["newList length lest than 25. Using all of cTable "]; newList = cTable[[endTerms ;;]]; ]; fitType = 2; fitFun[x_] = Switch[fitType, 1, Fit[Reverse@newList, {1, x}, x], 2, Fit[Reverse@newList, {1, x, x^2}, x], 3, Fit[Reverse@newList, {1, x, x^2, x^3}, x], 4, FindFormula[Reverse@newList, x, PerformanceGoal > "Quality"] ]; fitResults = fitFun[0] // N; Print[seriesIndex, ": ", N[fitResults, 10]]; lp2 = ListPlot[newList, PlotStyle > {PointSize[0.005], Red}]; plotA = Plot[fitFun[x], {x, 0, newList[[1, 1]]}, PlotStyle > {Dashed, Black}, AxesOrigin > {0, 0}]; xCoord = Plus @@ thePlotRange[[1]]/2; yCoord = (Plus @@ thePlotRange[[2]]/2); oneCyclePlot = Show[{cyclePlot, plotA, lp2}, PlotLabel > Style["1cycle series Root Test Results", 14, Bold, Black], PlotRange > {{0.001, 0.003}, {0.435, 0.455}}] oneCycleSummary = {seriesIndex, cycleType, Abs[theSingularList[[2]]] // N, fitResults // N, (Abs[fitResults  theSingularRegions[[2]]])/ theSingularRegions[[2]] 100 // N}
Each plot of Figure $3$ shows the scatter points of each $S$ set in blue, and the trailing set of points used in the fit shown as the red line of points. The red point on the vertical axis is $\lefts_2\right$. Its important to note that we determined one of the $2$cycle series continued onto a $1$cycle branch over $s_2$ and in fact continued until reaching a CLSP at $s_5$. However note that both series, as shown by the Root Test results are tending to the same $\lefts_2\right$ value as expected since the radius of convergence of a multicycled branch is the nearest CLSP.
The following is a summary of the Root Test results compared to the results determined by analytic continuation (AC) and the percent errors:
$$ \begin{array}{ccccc} \text{Series} & \text{Size} & \text{AC} & \text{Root Test} & \text{$\%$ error} \\ 2 & 1 & 0.438558 & 0.438787 & 0.0522727 \\ 1 & 2 & 0.438558 & 0.438886 & 0.0748519 \\ 3 & 2 & 0.438558 & 0.438886 & 0.0748519 \\ \end{array} $$
No comments:
Post a Comment