Note: This is a simple method of Newton polygon based on substitution which is useful for education purposes and generating series with a small number of terms. Background on the pologonal iteration step of the Newton method is covered in Section 6. Section 7 covers a more efficient method to generate series with many terms. Readers are encouraged to refer to Sections 6 and 7 for a practical and efficient method of Netwon Polygon.
The reader is asked to review the sections in sequential order as each succeeding section builds on a previous section.
In this part we use the principles described in the earlier sections to design a recursive function in Mathematica to compute all conjugate series of an algebraic function using the method of Newton Polygons. Our goal is to design a recursive function call as follows:
doPuiseux[currentFunction,functionIndex,totalIterations,normalExponent]
- currentFunction -- $F_n(z,w)$
- functionIndex -- $n$
- totalIterations -- total length of series
- normalExponent -- exponent used to find $w(z)$ for the original function $f(z,w)$
- theConjugateSeries -- global variable returning all conjugate series for function
theConjugateSeries={};
theFunction=z^4+2z^2w+(1+z^2)w^2+zw^3;
{normalizedFunction, lambda, mu} = convertToNormalForm[theFunction];
doPuiseux[normalizedFunction,1,10,lambda];
doPuiseux is to be called recursively for each conjugate series until we reach a polygon without a lower Newton leg or a lower leg consisting of a simple segment. In the later case, the recursion cycle stops and the remaining terms of the series are computed using substitution for a total series length of totalIterations. Our objective is to initially design the function only for cases in which we can solve the characteristic equation exactly. If we encounter the contrary condition, the function should terminate with a message that an exact solution cannot be found.
We wish to initially follow this flowchart:
And there are several conditions we wish to check for in order to make the algorithm robust.
- Residue terms in the series
- Plotting degenerate convex hulls
- Identifying non-exact roots to characteristic equation
Our objective is to work with exact arithmetic for the recursive stage in the algorithm and then numerically for the substitution stage. However, working with exact arithmetic can lead to zero-term large radical expressions which if not simplified, will give incorrect polygon results. For example, if we have $a_0(z)=k+z+z^2$ with $k$ being a very complex radical expression which is actually zero and if Mathematica does not recognize it is zero, when the convex hull is created for this function, the point $(0,0)$ will be used instead of the correct $(0,1)$. But simplifying every step in the analysis will significantly slow down execution time. However, Mathematica can very quickly determine the numerical value of $k$ but how do we determine if $k$ is really zero or just a number close to zero? At this point, we simply set a numerical threshold which is close to the numerical precision of the data we are working with. For example, if we are working with machine precision of $15$ decimal places and the numerical value of $k$ is $10^{-12}$, then we let $k$ equal exactly zero. This of course can lead to incorrect results if the value of $k$ is indeed close to zero but for now we accept that limitation. Therefore, during the analysis, we will check every coefficient this way and then if necessary, re-formulate the expressions with these "residue" terms removed.
The current Mathematica function use to plot convex hulls is ConvexHullMesh. However, this functions will not plot a set of co-linear points as is the case of degenerate polygons. But we would still like to have the means of plotting all the hulls of the analysis. We can work around this by first checking if the points are colinear. If they are, we simply detour around ConvexHullMesh and plot the points and the lines connecting them manually. We will create a function IsColinearQ[points] for this purpose.
Although we can design a method for working with non-exact roots to the characteristic equation, for now we wish to simply terminate the routine when this condition is encountered. Mathematica, if unable to solve for roots to a polynomial exactly, will represent the roots numerically as "root objects". The routine getNewSegment, will check for this condition and if found, end the routine. We use the following code for this:
Print["Solving for roots"];
mynetroots=x/.Map[Solve[#1==0,x]&,myeqnlist];
Print["Roots solved"];
If[Length[Cases[mynetroots,_Root,\[Infinity]]]>0,
{
Print[Style["Cannot find exact roots to characteristic equation",Red,Bold]];
Print["Char equ: ",myeqnlist];
Print["Roots: ",N[mynetroots]];
Abort[];
}
];
doPuiseux is a recursive routine and the reader is urged to review the sections describing the method. Two flags are passed to the routine: dFlag is used for diagnostic purposes. If True, then polygon, hull, characteristics equations, and other diagnostic information is reported during the analysis. Note in particular when doPuiseux is called recursively, we process the recursive function with removeFunctionResidues as described above. Often when the numerical results of a series is computed, the computations will leave small numerical residues. These are removed with a call to removeSeriesResidues.
This code contains color codes to make the code easier to follow. If imported into Mathematica, these codes must be removed.
doPuiseux[thef_,theK_,totalIter_,theLambda_,dFlag_,rootFlag_]:=
Module[{newtonPoints,hullPoints,gHull,gPoints,gnewtonLeg,gslopeLines,maxx,maxy,
slopeInterceptList,theCEquations,theRootTally,theclist,expList,minVals,currentSegment,
segmentnum,currentRootList,f,myroots,totalsegments,num,denom,num2,denom2,
theexp,mysegments,theSegmentData,theExponent,thePower,theDenom,tLength,p5,
rootnum,sum,totVal,myvals,elist,eindex,nindex,ctemp,atemp},
(* first check if initial function is reduceable *)
If[theK==1 && reduceQ[thef],
Print[Style["Starting function is reducible.",Red,Bold]];
Abort[];
];
theclist=CoefficientList[thef,w];
expList=Exponent[theclist,z,Min];
minVals=MapThread[Coefficient[#1,z,#2]&,{theclist,expList}];
(* get convex hull for this function *)
{newtonPoints,hullPoints,gHull,gPoints,gnewtonLeg,gslopeLines,maxx,maxy}=getPolygonSetup[thef,theK];
(* check if there is a lower newton leg *)
If[Length[newtonPoints]>1,
{
slopeInterceptList=getSlopeIntercepts[newtonPoints];
(* if in diagnostic mode, print convex hulls *)
If[dFlag,
{
Print["Lower Newton Leg points: ",newtonPoints];
Print["hullPoints: ",hullPoints];
(* check if this is a degenerate polygon *)
If[!isColinearQ[hullPoints],
{
p5=Graphics[Show[ConvexHullMesh[hullPoints]][[1]]/.
{ Directive[x_]:>Directive[{hexToRGB["#6495ed"],
Opacity[0.552],EdgeForm[{Black,Thickness[0.002]}]}]}];
Print[Show[{p5,gHull,gPoints,gnewtonLeg},
GridLines->{Table[i,{i,0,maxx}],Table[i,{i,0,maxy}]},AspectRatio->1,
PlotRange->{{0,maxx}, {0,maxy}},Axes->True,
PlotLabel->Style["Iteration " <>ToString[theK],16,Bold]]];
}
,
Print[Show[Graphics[{{PointSize[0.02],
Red,Point[hullPoints]},{Red,Line[hullPoints]}}],
Axes->True,PlotLabel->Style["Degenerate case detected",16]]];
];
Print["newtonPoints: ",newtonPoints];
];
(*get next segment in the lower leg*)
{slopeInterceptList,theCEquations,theRootTally}=getNewSegment[minVals,hullPoints,newtonPoints,theK];
If[rootFlag,
Print["root tally: ",theRootTally];
,
Print["root tally: ",N[theRootTally]];
];
theSegmentData=MapThread[{#1,#2,#3}&,{slopeInterceptList,theCEquations,theRootTally}];
If[dFlag,
Print["seg data: ",N[theSegmentData]];
];
mysegments=Range[Length[theSegmentData]];
If[dFlag,
Print["Total segments: ",mysegments];
];
(*for each segment, generate characteristic equation*)
Do[
{
If[dFlag,
Print["Segment Number: ",segmentnum];
];
Subscript[f, theK]=thef;
currentSegment=theSegmentData[[segmentnum]];
currentRootList=currentSegment[[3]];
(*assign lambda_k and beta_k from the segment data*)
{Subscript[\[Lambda], theK],Subscript[\[Beta], theK]}=currentSegment[[1]];
If[dFlag,
Print["lambda ",theK," is: ",Subscript[\[Lambda], theK]," and beta ",
theK," is: ",Subscript[\[Beta], theK]];
];
(*for each root of characteristic equation, generate a power series*)
For[rootnum=1,rootnum<=Length[currentRootList],rootnum++,
Subscript[c, theK]=currentRootList[[rootnum,1]];
If[dFlag,
Print["thek: ",theK];
Print["c",theK," is ", N[Subscript[c, theK]]];
];
Subscript[f, theK+1]=Expand[1/z^Subscript[\[Beta],
theK] Subscript[f, theK]/.w->z^Subscript[\[Lambda],theK] (Subscript[c, theK]+w)];
Print[polyForm[N[Subscript[f, theK+1]],w]];
(*check if multiple root then enter next recursive cycle*)
If[currentRootList[[rootnum,2]]!=1 || Subscript[\[Lambda], theK]==0,
{
Print["Entering recursion level ",theK+1," with f: ",
N[polyForm[Subscript[f, theK+1],w]]];
doPuiseux[removeFunctionResidues[Subscript[f,theK+1]],
theK+1,totalIter,theLambda,dFlag,rootFlag];
}
,
(*do regular substitution here*)
{
(*first get denominator of last lambda*)
theDenom=Denominator[Subscript[\[Lambda], theK]];
Print["and denum: ",theDenom];
sum=0;
(*first sum up the values from the recursive steps*)
For[totVal=1,totVal<=theK,totVal++,
sum+=N[Subscript[c, totVal]] z^Sum[Subscript[\[Lambda], i],{i,1,totVal}];
];
If[dFlag,
Print["initial sum: ",sum];
];
(*next, initialize the subscripts for the coefficients to undefined*)
Quiet[Table[Subscript[a, n]=.,{n,1,20}]];
(*now form the place-holders for the sum a_(theK+n) z^(n/theDenom)*)
myvals=Plus @@ Table[Subscript[a, theK+n] z^(n/theDenom),{n,1,totalIter-theK}];
If[dFlag,
Print["myvals: ",myvals];
];
(*substitute myvals into the function*)
Subscript[f, theK+1]=Expand[N[Subscript[f, theK+1]/.w->myvals]];
(*get exponent list *)
elist=Exponent[Subscript[f, theK+1],z,List]//Sort;
If[dFlag,
Print["the exponent list: ",elist];
];
eindex=1;
nindex=theK+1;
(*solve for each of the unknown coefficients a_n*)
While[nindex<=totalIter,
{
ctemp=Coefficient[Subscript[f, theK+1],z,elist[[eindex]]];
atemp=Coefficient[ctemp,Subscript[a, nindex]];
While[atemp==0,
{
eindex++;
ctemp=Coefficient[N[Subscript[f, theK+1]],z,elist[[eindex]]];
atemp=Coefficient[ctemp,Subscript[a, nindex]];
}
];
Subscript[a, nindex]=N[Subscript[a, nindex]/.NSolve[ctemp==0,Subscript[a, nindex]]];
Subscript[\[Lambda], nindex]=1/theDenom;
sum+=Subscript[a, nindex] z^Sum[Subscript[\[Lambda], i],{i,1,nindex}];
If[dFlag,
Print["nval: ",nindex, " aindex: ",Subscript[a, nindex]];
];
nindex++;
eindex++;
}
];
(*now print the series*)
Print[Style["Series ",Black, Bold],
Style[Length[theFunctionSeries]+1,Black,Bold],
Style[" is: ",Black,Bold],Style[Expand[removeSeriesResidues[sum]/(z^theLambda)],Black,Bold]];
(*append the series to the FunctionSeries list*)
theFunctionSeries=Append[theFunctionSeries,Expand[removeSeriesResidues[sum]/(z^theLambda)]];
}
];
];
},
{segmentnum,mysegments}
];
}
,
(*at this point, we have no lower leg*)
{
Print["length of newtonPoints is less than two. Series is terminating. "];
Print["newtonPoints: ",newtonPoints];
sum=0;
(* sum the finite series *)
For[totVal=1,totVal<=theK-1,totVal++,
sum+=N[Subscript[c, totVal]] z^Sum[Subscript[\[Lambda], i],{i,1,totVal}];
];
Print[Style["Series ",Black, Bold],
Style[Length[theFunctionSeries]+1,Black,Bold],
Style[" is: ",Black,Bold],Style[Expand[removeSeriesResidues[sum]/(z^theLambda)],Black,Bold]];
theFunctionSeries=Append[theFunctionSeries,Expand[removeSeriesResidues[sum]/(z^theLambda)]];
}
];
];
We would for example call the routine as follows:
theFunction=(z+z^2)+(z+z^2)w+(z+z^3)w^2+(z) w^3+(1+z)w^4;
theFunctionSeries={};
totalIterations=10;
Quiet[Table[Subscript[a, n]=.,{n,1,totalIterations}]];
{normalizedFunction,lambda,mu}=convertToNormalForm[theFunction];
doPuiseux[normalizedFunction,1,totalIterations,lambda,False,False]
No comments:
Post a Comment