Mathematica Code

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)} $$
  1. Initialization functions


    1. {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];
      
      
    2. {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 plot
      

      Function 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}
      ]
      
      
    3. 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
      ];
      
      
      
    4. 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
      ];
      
      
      
      
    5. 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]; 
      
      
      
    6. {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}
      
      
      ];
      
      
  2. 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.


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

Blog Archive