Refer to the pages above for more information about these functions. This page is currently under construction. When completed it will contain the complete code for implementing the Newton Polygon algorithm in Mathematica.
$$ \newcommand{\bint}{\displaystyle{\int\hspace{-10.4pt}\Large\mathit{8}}} \newcommand{\res}{\displaystyle{\text{Res}}} \newcommand{\wvalx}{\underbrace{z^{\lambda_4}(c_4+w_5)}_{w_4}} \newcommand{wvalxx}{\underbrace{z^{\lambda_3}(c_3+\wvalx)}_{w_3}} \newcommand{wvalxxx}{\underbrace{z^{\lambda_2}\{c_2+\wvalxx\}}_{w_2}} \newcommand{wvalxxxx}{z^{\lambda_1}\big(c_1+\wvalxxx\big)} $$
-
Initialization functions
-
{normalizedFunction,lambda,mu}=getNormalizedFunction[function]This routine normalizes a function. For example:
theFunction=(1+z)+(z+z^2)w+z^3 w^2; {normalizedFunction,lambda,mu}=getNormalizedFunction[theFunction];
returns{(z+z^2)+(1+z)w+w^2,2,1}
Function code
convertToNormalForm[p_]:=Module[{mymu,temp,flag,mylambda,thepoly,degree ,currentmu,mychecks,mynegs,theclist,theOrders}, theclist=CoefficientList[p,w]; degree=Length[theclist]-1; mylambda=0; theOrders=Table[getTheOrder[theclist[[n+1]]],{n,0,degree}]; If[theOrders[[degree+1]]!=0,{mymu[n_,\[Lambda]_]:=n \[Lambda]-theOrders[[degree+1]]; flag=False; While[flag==False,mylambda++; currentmu=mymu[degree,mylambda]; If[currentmu>0,{mychecks=Table[currentmu+theOrders[[i]]-(i-1) mylambda,{i,2,degree}]; mynegs=Select[mychecks,#<0&]; If[Length[mynegs]>0,flag=False;,flag=True;];},flag=False;];]; thepoly=Expand[z^currentmu (p/.w->w/(z^mylambda))];},{thepoly=p;}]; {thepoly,mylambda,currentmu}]; getTheOrder[p_]:=Module[{clist,theorder,j},clist=CoefficientList[p,z]; theorder=\[Infinity]; If[Length[clist]!=0,{j=0; While[clist[[j+1]]==0&&(j+1)<=Length[clist],{j++;}]; theorder=j;},theorder=\[Infinity];]; theorder]; -
{newtonPoints, hullPoints, gHull, gPoints, gnewtonLeg, gslopeLines, maxx, maxy} = getPolygonSetup[normalizedFunction, type]; (and associated functions)Returns the Newton polygon for a normalized function. If type=1 then segments with zero slope are included with the returned data otherwise only legs with negative slopes are returned. Example:
theFunction=(1+z)+(z+z^2)w+z^3 w^2; {normalizedFunction,lambda,mu}=getNormalizedFunction[theFunction]; {newtonPoints, hullPoints, gHull,gPoints, gnewtonLeg, gslopeLines, maxx, maxy} = getPolygonSetup[normalizedFunction, 1];
returns:newtonPoints={point1,point2,...,pointn} (* points on the lower Newton leg *) hullPoints={hull1,hull2,...,hullm} (* all points making up the convex hull *) gHull=graphics polygon of hull gPoints=graphics points of data set gnewtonLeg=graphics of Newton leg gslopeLines=graphics of slopes of segments maxx=maximum x for plot maxy=maximum y for plotFunction code
Needs["ComputationalGeometry`"]; hexToRGB=RGBColor@@(IntegerDigits[ToExpression@StringReplace[#,"#"->"16^^"],256,3]/255.)& polyForm[poly_,var_]:=Module[{coeffs=CoefficientRules[poly,var]//Sort},Interpretation[Row[Table[Row[{"(",coeff[[2]],")",var^coeff[[1,1]]/.1->""}],{coeff,coeffs}],"+"],poly]] myyLine[x_,x1_,y1_,x2_,y2_]:=(y2-y1)/(x2-x1) (x-x1)+y1; lowerConvexHull[ls_]:= Module[{ v=ConvexHull[ls],ks,min,max,left,right,pp},ks=ls[[v]]; min=Min[ks[[All,1]]]; pp=Position[ks,{min,_}]; left=Last[pp][[1]]; Append[Take[ks,left-Length[ks]-1],First[ks]]]; getPolygonSetup[fn_,thek_]:=Module[{mypoints,lcp,lcpsave,slope,mylines,newtonleg,sg,convexhull,pgp,maxx,maxy,wdegree}, wdegree=Exponent[fn,w]; mypoints=CoefficientRules[fn ,w]/.Rule[u_,v_]:> {First[u],Exponent[v,z,Min[#&]]}; lcp=lowerConvexHull[mypoints]; lcpsave=lcp; For[n=1,n<=Length[lcp]-1,n++, slope=(lcp[[n+1,2]]-lcp[[n,2]])/(lcp[[n+1,1]]-lcp[[n,1]]); If[thek==1, { If[slope>0, lcpsave=DeleteCases[lcpsave,lcp[[n+1]]]; ]; } , { If[slope>=0, lcpsave=DeleteCases[lcpsave,lcp[[n+1]]]; ]; } ]; ]; mylines=Table[Plot[myyLine[x,lcpsave[[n,1]],lcpsave[[n,2]],lcpsave[[n+1,1]],lcpsave[[n+1,2]]],{x,0,wdegree},PlotStyle->{Dashed,Blue}],{n,1,Length[lcpsave]-1}]; newtonleg=Graphics[{Red,Thickness[0.007],Line[lcpsave]},Axes->True]; (* for fully-ramified branch in normal form, newton polygon is from first to last point *) sg=Graphics[{Red,PointSize[0.02],Point @@{mypoints}}]; convexhull=ConvexHull[mypoints]; maxy=Max[(#[[2]]& /@ mypoints)]; maxx=Max[(#[[1]]& /@ mypoints)]; If[Length[mypoints]>2, (* cannot draw convex hull when less than 3 points *) pgp=PlanarGraphPlot[mypoints,convexhull,Axes->True,PlotLabel->"Newton Polygon",PlotRange->{{0,15},{0,10}}]; , pgp=Graphics[Null]; ]; {lcpsave,mypoints,pgp,sg,newtonleg,mylines,maxx,maxy} ] -
theSlopeInterceptList=getSlopeIntercepts[newtonPoints];Returns the slope and intercept for each segment on the lower Newton leg. For example:
theFunction=(1+z)+(z+z^2)w+z^3 w^2; {normalizedFunction,lambda,mu}=getNormalizedFunction[theFunction]; {newtonPoints, hullPoints, gHull,gPoints, gnewtonLeg, gslopeLines, maxx, maxy} = getPolygonSetup[normalizedFunction, 1]; theSlopeIntercepts=getSlopeIntercepts[newtonPoints];
returns:{{slope1,intercept1},{slope2,intercept2},...,{slopen,interceptn}}Function code
getSlopeIntercepts[lowerLeg_]:=Module[{lnp,myeqns,myslopes,myinter,mypointlist,slopeinterceptlist}, lnp=lowerLeg; myeqns=Table[myyLine[x,lnp[[n,1]],lnp[[n,2]],lnp[[n+1,1]],lnp[[n+1,2]]],{n,1,Length[lnp]-1}]; myslopes=Coefficient[myeqns,x]; myinter=myeqns/.x->0; mypointlist=MapThread[{#1,#2}&,{-myslopes,myinter}]; slopeinterceptlist=DeleteDuplicates[mypointlist]; slopeinterceptlist ]; -
coefficientMatrix=getCoefficientMatrix[function];Returns a matrix of exponents and coefficients. The exponents on z may be fractional. Each row of the matrix has an exponent list and coefficient list. The rows are the coefficients of w^n starting with n=0. For example, suppose we are given: $$ \begin{aligned} \text{theFunction}&=(1+2c^3+c^6+c^2 z^{1/3}+c^3 z^5) \\ &+(6 c^2+6 c^5+2 c z^{1/3}+3 c^2 z^5)w \\ &+(6 c+15 c^4+z^{1/3}+3 c z^5) w^2 \\ &+(2+20 c^3+z^5) w^3 \\ &+ 15 c^2 w^4 \\ &+ 6c w^5 \\ &+ w^6 \end{aligned} $$
then:coefficientMatrix=getCoefficientMatrix[theFunction]
returns:{{{0,1/3,5},{1+2 c^3+c^6,c^2,c^3}}, {{0,1/3,5},{6 c^2+6 c^5,2 c,3 c^2}}, {{0,1/3,5},{6 c+15 c^4,1,3 c}}, {{0,5},{2+20 c^3,1}}, {{0},{15 c^2}}, {{0},{6 c}}, {{0},{1}}}Function code
getCoefficientMatrix[thef_]:=Module[{plist,coefficientMatrix,exponents,zlist,array,temp}, plist=CoefficientList[thef,w]; coefficientMatrix={}; For[i=1,i<=Length[plist],i++, exponents=Exponent[plist[[i]],z,List]; zlist=Coefficient[plist[[i]],z,exponents]; coefficientMatrix=Append[coefficientMatrix,{exponents,zlist}]; ]; coefficientMatrix ]; -
characteristicEquations=getCharacteristicEquations[function,index];Returns a list of characteristic equations for a Newton polygon. If index=1, returns negative and zero slope segments. Otherwise, returns only negative slope segments. For example, suppose: $$ F(z,w)=(z^4)+(2 z^2+z^4)w+(1+z^2+z^3)w^2+(z)w^3+(1/4-z/2)w^4+(-(1/2))w^5 $$
then:{lcpsave,myPoints,pgp,sg,newtonLeg,mylines,maxx,maxy}=getPolygonSetup[normalizedFunction,1] slopeInterceptList=getSlopeIntercepts[lcpsave] theftemp=theFunction; theclist=CoefficientList[theftemp,w]; expList=Exponent[theclist,z,Min]; minVals=MapThread[Coefficient[#1,z,#2]&,{theclist,expList}]; theEquations=getCharacteristicEquations[minVals,myPoints,lcpsave,1]
returns:{{1+2c_1+c_1^2=0, c_1^2+1/2 c_1^4-1/2 c_1^5=0}Function code
getCharacteristicEquations[minvals_, mypoints_, lnp_, thek_] := Module[{myeqns, myslopes, myinter, mypointlist, slopeinterceptlist, totalsegmentsize, pointlist, myeqnlist, myeqn}, myeqns = Table[myyLine[x, lnp[[n,1]], lnp[[n,2]], lnp[[n + 1,1]], lnp[[n + 1,2]]], {n, 1, Length[lnp] - 1}]; myslopes = Coefficient[myeqns, x]; myinter = myeqns /. x -> 0; mypointlist = MapThread[{#1, #2} & , {-myslopes, myinter}]; slopeinterceptlist = DeleteDuplicates[mypointlist]; Print["slopeinterceptlist: ", slopeinterceptlist]; totalsegmentsize = Length[slopeinterceptlist]; pointlist = {}; For[i = 1, i <= Length[slopeinterceptlist], i++, pointlist = Append[pointlist, {}]; eqn[x_] := (-slopeinterceptlist[[i,1]])*x + slopeinterceptlist[[i,2]]; For[j = 1, j <= Length[mypoints], j++, If[eqn[mypoints[[j,1]]] == mypoints[[j,2]], pointlist[[i]] = Append[pointlist[[i]], mypoints[[j]]]; ]; ]; ]; myeqnlist = {}; For[i = 1, i <= Length[pointlist], i++, myeqn = Sum[minvals[[1 + pointlist[[i,j,1]]]]*x^pointlist[[i,j,1]], {j, 1, Length[pointlist[[i]]]}] == 0; myeqnlist = Append[myeqnlist, myeqn]; ]; myeqnlist]; -
{slopeIntercept,characteristicEquation,rootTally}=
getNewSegment[minVals, myPoints, lcpsave, theK];For a Newton polygon, returns for each segment on the lower Newton leg, the slope and intercept of the segment, the associated characteristic equation, and the roots to the characteristic equation along with their multiplicities. For example, suppose: $$ F(z,w)=(z^4)+(2 z^2+z^4)w+(1+z^2+z^3)w^2+(z)w^3+(1/4-z/2)w^4+(-(1/2))w^5 $$
then:theFunction=w^2+2z^2 w+z^4+z^2 w^2+z w^3+1/4 w^4+z^4 w+z^3 w^2-1/2 z w^4-1/2 w^5; theclist=CoefficientList[theFunction,w]; expList=Exponent[theclist,z,Min]; minVals=MapThread[Coefficient[#1,z,#2]&,{theclist,expList}]; {lcpsave,myPoints,pgp,sg,newtonLeg,mylines,maxx,maxy}=getPolygonSetup[theftemp,theK]; {slopeInterceptList, charEquations, theRootTally} = getNewSegment[minVals, myPoints, lcpsave, 1]; theSegmentData=MapThread[{#1,#2,#3}&,{slopeInterceptList,theCEquations,theRootTally}];
returns as theSegmentData:{{{2,4}, 1+2 x+x^2==0, {{-1,2}}}, {{0,0}, x^2+x^4/4-x^5/2==0, {{0,2}, {1/6 (1+(217-12 Sqrt[327])^(1/3)+(217+12 Sqrt[327])^(1/3)),1}, {1/6-1/12 (1+I Sqrt[3]) (217-12 Sqrt[327])^(1/3)-1/12 (1-I Sqrt[3]) (217+12 Sqrt[327])^(1/3),1}, {1/6-1/12 (1-I Sqrt[3]) (217-12 Sqrt[327])^(1/3)-1/12 (1+I Sqrt[3]) (217+12 Sqrt[327])^(1/3),1}}}}Where: $\{2,4\}$ is the slope and intercept of the first segment, $1+2x+x^2=0$ is the associated characteristic equation, and $\{\{-1,2\}\}$ is the root tally, that is, $-1$ with multiplicity 2. For the horizontal segment we next have $\{0,0\}$ as the slope and intercept, $x^2+x^4/4-x^5/2=0$ as its characteristic equation, and $\{0,2\}$ as a root of zero with multiplicty two and the three remaining roots with multiplicity one
Function code
getNewSegment[minvals_,mypoints_,lnp_,thek_]:=Module[{myeqns,myslopes,myinter,mypointlist,slopeinterceptlist,totalsegmentsize,pointlist,myeqnlist,myeqn,mynetroots,roottally,workingPrecision,rootmatrix,rootaccuracy,rootprecision,ra,eqnlistprecision,workingprecision,eqnlistPrecision,eqnlistAccuracy,minvalaccuracy,myeqnstatus,temp,rootstatuslist,rootstatus,rootprecisionarray,rootsummary}, myeqns=Table[myyLine[x,lnp[[n,1]],lnp[[n,2]],lnp[[n+1,1]],lnp[[n+1,2]]],{n,1,Length[lnp]-1}]; myslopes=Coefficient[myeqns,x]; myinter=myeqns/.x->0; mypointlist=MapThread[{#1,#2}&,{-myslopes,myinter}]; slopeinterceptlist=DeleteDuplicates[mypointlist]; Print["slopeinterceptlist: ",slopeinterceptlist]; totalsegmentsize=Length[slopeinterceptlist]; (* now have slope-intercept list for each lower segment. Now check each point to see what points are on those lines *) pointlist={}; For[i=1,i<=Length[slopeinterceptlist],i++, pointlist=Append[pointlist,{}]; eqn[x_]:=-slopeinterceptlist[[i,1]] x+slopeinterceptlist[[i,2]]; For[j=1,j<=Length[mypoints],j++, If[eqn[mypoints[[j,1]]]==mypoints[[j,2]], pointlist[[i]]=Append[pointlist[[i]],mypoints[[j]]]; ]; ]; ]; myeqnlist={}; For[i=1,i<=Length[pointlist],i++, myeqn=\!\( \*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(Length[pointlist[\([\)\(i\)\(]\)]]\)]\(minvals[\([\)\(1 + pointlist[\([\)\(i, j, 1\)\(]\)]\)\(]\)]\ x^\((pointlist[\([\)\(i, j, 1\)\(]\)])\)\)\); If[Resultant[myeqn,D[myeqn,x],x]==0, Print[Style["Multiple roots detected",Blue]]; ]; myeqnlist=Append[myeqnlist,myeqn]; ]; myeqnstatus={}; 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[]; } ]; roottally=DeleteCases[#,{0,_}]& /@ (Map[Tally[#1]&,mynetroots]); {slopeinterceptlist,myeqnlist,roottally} ];
-
-
doPuiseux[thef_, theK_, totalIter_, theLambda_, dFlag_,rootFlag_]
Computes the fractional power series for the function $w(z)$ given implicitly as $$ f(z,w)=a_0(z)+a_1(z)w+a_2(z)w^2+\cdots+a_n(z)w^n=0 $$ Please refer to the section on doPuiseux for a description of this function.
-
getDigits[{vals}_, accuracy_]
Returns an integer representation of a floating-point number with accuracy number of digits.
getDigits[vals_, n_] := Module[{net, theDigits, i, theChoppedVals, theComponents, myList, listAccuracy}, If[Length[vals] == 0, myList = {vals}; , myList = vals; ]; listAccuracy = Accuracy[myList]; myList = SetAccuracy[myList, listAccuracy - 2]; net = Length[myList]; theComponents = {Re[#], Im[#]} & /@ myList; theChoppedVals = Chop[theComponents, 10^(-n)]; theDigits = RealDigits[theChoppedVals, 10, n]; For[i = 1, i <= net, i++, If[theChoppedVals[[i, 1]] < 0, theDigits[[i, 1]] = Prepend[theDigits[[i, 1]], -1]; , theDigits[[i, 1]] = Prepend[theDigits[[i, 1]], 1]; ]; If[theChoppedVals[[i, 2]] < 0, theDigits[[i, 2]] = Prepend[theDigits[[i, 2]], -1]; , theDigits[[i, 2]] = Prepend[theDigits[[i, 2]], 1]; ]; ]; theDigits ];
References:
1 Algebraic Curves, Walker, R.J.2 All Algebraic Functions can be Computed Fast, H. T. Kung and J. F. Traub
No comments:
Post a Comment