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 built-in function $\texttt{AsymptoticSolve}$ can compute the expansions although there are some issues. In order to determine the radius of convergence of a power expansion, impinging singular points (ISPs) are identified for each branch sheet over a base singular point $s_b$ by the method of analytic continuation: the sheets are numerically-integrated across successively distant singular points until each branch sheet analytically continues onto either a multi-valued branch sheet or a single valued sheet with a pole. The singular point where this occurs is the ISP for the branch sheet. Different branch sheets of a multi-cycled branch will in general have different ISPs. The ISP nearest to the expansion center is the convergence-limiting singular point (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 analytically-continuous 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 analytically-continuous across successively distant singular points over analytic $1$-cycle branches until the first ISP is encountered idenfified by $s_c$. This singular point becomes the CLSP and establishes the radius of convergence, $R$ of the branch expansions given by $R=\left|s_b-s_c\right|$
- 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 over-written by the black dots.
Mathematica code 1
$$ \begin{array}{c} \\ \text{Table 1: Singular points for (2)} \end{array} \\ \begin{array}{cc} 1 & 0 \\ 2 & -0.3582-0.2530 i \\ 3 & -0.3582+0.2530 i \\ 4 & 0.1961-0.5259 i \\ 5 & 0.1961+0.5259 i \\ 6 & \color{red}{0.6823} \\ 7 & 0.7492 \\ 8 & -0.1440-0.9401 i \\ 9 & -0.1440+0.9401 i \\ 10 & \color{red}{-0.3412-1.1615 i} \\ 11 & \color{red}{-0.3412+1.1615 i} \\ 12 & -0.843-1.560 i \\ 13 & -0.843+1.560 i \\ \end{array} $$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
The series are next sorted relative to each singular point reference point $s_n+1/3r$. In this example, for low-degree functions and with a rational center, $\texttt{AsymptoticSolve}$ can solve for the exact coefficients of the Puiseux expansions. The first $8$ terms of the series computed in code block $2$ are $$ \begin{align*} P_1(z)&=-1. z^{1/2}-1. z^{3/2}+0.5 z^2-1. z^{5/2}+1.5 z^3-1.625 z^{7/2}+z^4-2.625 z^{9/2}\\ P_2(z)&=z^2-1. z^3+4. z^7-9. z^8+9. z^9-7. z^{10}-12. z^{11}+91. z^{12}\\ P_3(z)&=z^{1/2}+z^{3/2}+0.5 z^2+z^{5/2}+1.5 z^3+1.625 z^{7/2}+z^4+2.625 z^{9/2} \end{align*} $$ where numerical values of the coefficients are given for brevity. By inspection, we see a $2$-cycle conjugate class and a single-cycle series. Since the lowest non-zero exponent of the $2$-cycle is $1/2$, this is a $V_2$ branch and the single-cycle is a Taylor series for a $T$ branch. The sorted cycle sizes in this example is $\{2,1,2\}$ representing the cycle size of each Puiseux series in sorted order. Note: if there were more than $k$ series of $k$-cycles, additional notation would be needed to identify which series belongs to which conjugate class. See reference above.(* 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 // MatrixFormHowever in the case of a finite-precision 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 analytically-continuous 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.3852-0.2530i,-0.3852+0.2530i\}$ and our objective is to create an analyticallly-continuous 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(x-x_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*} (x-x_b)^2+(y-y_b)^2&=r_b^2 \\ (x-x_n)^2+(y-y_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(1-t)+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.3582-0.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(1-t)+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(A-s_b); \quad i=1,2,3 \end{equation} $$ with $z(t)=\text{pointA}(1-t)+\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.3582-0.2530i$ up to degree $512$ and sort:
Mathematica code 3
The first $5$ terms of the (sorted) series at $s_2$ are: $$ \begin{align*} Q_1(z)&=(0.030+0.354 I)+(0.073+0.537 I) z^{1/2}-(0.629+0.264 I) z-(0.560+0.389 I) z^{3/2}+(0.082-0.501 I) z^2 \\ Q_2(z)&=(0.030+0.354 I)-(0.073+0.537 I) z^{1/2}-(0.629+0.264 I) z+(0.560+0.389 I) z^{3/2}+(0.082-0.501 I) z^2\\ Q_3(z)&=(0.095-0.476 I)+(0.175+0.312 I) z+(0.210-0.002 I) z^2+(0.543-0.089 I) z^3+(1.034-0.582 I) z^4 \end{align*} $$ so that the branch sequence over this singular point is $\{2,2,1\}$. So that we have the following branching over $s_b$ and $s_2$: $$ \begin{array}{ccc} \text{Series Number} & s_b & s_n \\ 1 & 2 & 2 \\ 2 & 1 & 2 \\ 3 & 2 & 1 \\ \end{array} $$ And consider the necessary and sufficient conditions for branch continuations described in the paper: An $n$-cycle branch can continue across a singular point $s_n$ only if $s_n$ has $n$ analytic single-cycle branches, that is, $n$ single-cycle branches which are not poles. So that since there are only one $1$-cycle branch over $s_n$, the $2$-cycle branch over $s_b$ cannot continue across $s_b$. The only possible continuation across $s_n$ is if the $1$-cycle branch at $s_b$ continues onto the $1$-cycle branch over $s_n$. We now solve the IVPs above to determine this: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 // MatrixFormMathematica code 4
So that we integrated over the branch sheets $\{2,1,2\}$ of $s_b$ onto point $D$ in Figure $2$ of the branch sheets of $s_2$ returning the values $\{-0.0471632 + 0.43161 I, 0.0218879 + 0.105942 I, 0.0905285 - 0.416981 I\}$. We next compute the value of each $\text{$Q_i$(pointD-thesingularList[[2]])}$ and obtain $\{-0.0471632 + 0.43161 I, 0.0218879 + 0.105942 I, 0.0905285 - 0.416981 I\}$ and we summarize the results: $$ \begin{array}{cccc} \text{$s_b$ branch type} & \text{$w_i$(1)} & \text{$Q_i$(pointD-$s_2$)}&\text{$s_2$ branch type} \\ 2& -0.04716+0.4316i & -0.04716+0.4316i & 2\\ 1 & 0.02188+0.1059i & 0.02188+0.1059i & 2 \\ 2 & 0.0905-0.4169i & 0.09052-0.4169i & 1 \end{array} $$ Consider the first line in the table: We analytically continued a $2$-cycle sheet over $s_b$ onto the $s_2$ domain at point $D$ with value $-0.04716+0.4316i$. And we expanded the branching over $s_2$ into three series with cycle orders $\{2,2,1\}$. The first series over $s_2$ is a sheet of a $2$-cycle and $\text{$Q_1$(pointD-$s_2$)}=-0.04716+0.4316i$. This means the $2$-cycle sheet over $s_b$ continued onto a $2$-cycle sheet of $s_2$. And thus the $2$-cycle branch over $s_b$ is not continuous over $s_2$. Now note the second line in the table: the $1$-cycle sheet over $s_b$ continued onto the other $2$-cycle sheet of $s_2$ so is also not continuous over $s_2$. And finally note in the third line of the table that even though the second $2$-cycle of $s_b$ continued onto a $1$-cycle over $s_2$, the radius of convergence of the branch is the nearest impinging singular point of the branch sheets. We have therefore reached the CLSP for the series over $s_b$ and can conclude the (exact) radius of convergence for the three series over $s_b$ is $\left|s_2\right|$.zTrace[t_]:=pointA(1-t)+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->(nextNDZ-theSingularList[[2]])//NThis 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 a branch sheet reaches the CLSP.
-
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=\left|s_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["1-cycle 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 $\left|s_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 $\left|s_2\right|$ value as expected since the radius of convergence of a multi-cycled 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