Gauss Topology  October 2015

By  Barry H Dayton, www.barryhdayton.space

Available Procedures

Important Note :  Although many procedures ask for variables X, the plotting procedures only work for X = {x, y}.

GaussCurve[p]  Finds the Gauss Curve associated with the complex polynomial  p(z)

QPlot[f, r, X]   Gives a quick plot of the curve without decomposition, showing the region f(x,y)>0 inside the circle of radius r about the origin.  No decomposition is made.

QData[f, r, X]  Gives the raw data assoiated with QPlot, a list of lists of points on the curve, organized in no particular pattern.

GPlot[f ,r, X]  Plots the Gauss Decomposition of f inside circle of radius r

GData[f, r, X, v]  Data associated with GPlot, each component is a list of points.

GListPlot[f, L]  Plots the user list or lists of points L on background of region f(x,y)>0 inside a circular region big enough only to include points of L.  An optional 3rd argument is a list of color data.

RefineComponent[f, L, X]   Refines a component, approximately doubling the number of points in the list.  Works on only one component at a time.  Warning: L must match f in direction.

FindComp[f, p, L,  X]  Attempts to find the component of the decomposition of f given by L and add p to this component in correct spot.

See also the notebook  GaussTopExamples.nb  for examples and further options.

Code

GaussCurve[h_]:=ComplexExpand[Re[h/.{z->x+I*y}]]

X={x,y}

tangentrealpoints[f_,kk_,X_,Q_:True]:=Module[{J, sol,rp, p,q,k},
If[(kk/.Thread[X->RandomReal[{-1,1},2]])==0, k=RandomReal[{-1,1},2].X, k=kk];
If[Q,Print["k= ",k]]; 
J=Det[D[{f,k},{X}]];
sol=NSolve[{f,J}] ;
rp=Cases[X/.sol,{_Real,_Real}];
If[rp=={}&& Q,Print["trp: no point found"];Return[]];
rp
]

line[p_,q_,X_]:=Module[{e1,e2,a,b,c,M},
e1=Switch[Length[p],2,Append[p,1],3,p];
e2=Switch[Length[q],2,Append[q,1],3,q];
M={e1,e2,RandomReal[{-1,1},3]};
({a,b,c}/.(Solve[M.{a,b,c}=={0,0,1}][[1]])).Append[X,1]
]

tline[f_,p_,X_]:=line[p,{D[f,X[[2]]],-D[f,X[[1]]],0}/.Thread[X->p],X]

realsingularpoints[f_,X_,eps_:.02]:=
Flatten[Take[Select[Tally[Cases[X/.NSolve[{f,{D[f,X[[1]]],D[f,X[[2]]]}.RandomReal[{-1,1},2]}],{_Real,_Real}],Norm[#1-#2]<eps&],#[[2]]>1&],All,1],1]

(* Total degree of multivariate polynomial  *)
tdeg[f_,X_]:= With[{T= Thread[X->Table[t_,{i,Length[X]}]]}, 
Max[Exponent[MonomialList[f]/.T,t_]]
]

(* problem -- seems only to work if X={x,y}  *)
realinflectionpoints[f_,X_]:=Module[{H,fh,d,tt,rp,L,X1,sol,p},
d=tdeg[f,X];
fh=Expand[tt^d*(f/.Thread[X->1/tt*X])];
X1=Append[X,tt];
H=Det[D[fh,{X1,2}]];
If[H==0,Return[{}]];
rp=Cases[X/.NSolve[{fh,H,tt-1}],{_Real,_Real}];
If[Length[rp]>0,L= rp,Return[{}]];
rp=Reap[Do[
If[Norm[{D[f,X[[1]]],D[f,X[[2]]]}/.Thread[X->L[[i]]]]> .01, Sow[L[[i]]]], {i,Length[L]}]][[2]];
If[Length[rp]>0,rp[[1]],{}]
];

infiniterealpoints[f_,X_]:=Module[{F,pp,sol,RR,iptr,hr, hp, ml,rp,k, degf,ips},
ml=MonomialList[f];
degf=tdeg[f,X];
hp=Total[Select[ml,tdeg[#,X]== degf&]];
k=RandomReal[{-1,1},2].X-1;
sol=NSolve[{hp,k}];
rp=Cases[X/.sol,{_Real,_Real}];
If[rp=={},Return[{}],
ips=rp.{{1,0,0},{0,1,0}}];
Normalize/@DeleteDuplicates[ips]
];

GFrame[f_,p_,X_]:=Module[{tl,nl,tv,nv,u},
 tv={D[f,X[[2]]],-D[f,X[[1]]]}/.Thread[X->p];
If[Norm[tv]>.001,tv=Normalize[tv],Return[{0,0}]];
nv=Cross[tv];
u=p+.1*nv;
If[(f/.Thread[X->u])<0,nv=-nv];
tv=Cross[nv];
{tv,nv}]

closestrealpoint[f_,p_,X_]:=Module[{sol},
sol=Cases[X/.NSolve[{{D[f,X[[2]]],-D[f,X[[1]]]}.{X[[1]]-p[[1]],X[[2]]-p[[2]]},f}],{_Real,_Real}];
First[SortBy[sol,Norm[#-p]&]]
]

almostMemberQ[S_,p_,ep_]:=Length[Select[S,Norm[#-p]<2*ep&]]>0

npos[L_,a_,eps_]:=Module[{len,flag},
len=Length[L];
flag=0;
Do[If[Norm[L[[i]]-a]<eps,flag=i;Break[]],{i,len}];
If[flag==0,Return[{}],Return[{{flag}}]]
]

GCP[f_,r_,v_,X_]:=Module[{p,q,sol, bp,gbp,gcp,cp,tp,tpx,tpy,tal,ip,tv,gep,del,stv,sp,w,nsp,sol2},
del=  0.2;
bp=Cases[X/.NSolve[{f,X[[1]]^2+X[[2]]^2-r^2}],{_Real,_Real}];
gbp={};
gep={};
tp={};
sp=Select[realsingularpoints[f,X,.022],Norm[#]<r+.001&];
Do[p=bp[[i]];
tv=First[GFrame[f,p,X]];
If[almostMemberQ[sp,p,.022],Print["Singular Boundary Point, increase r"];Abort[]];
 If[Norm[p+tv]>Norm[p], gep=Append[gep,p], gbp=Append[gbp,p]],
{i,Length[bp]}];
tpx=Select[tangentrealpoints[f,v[[1]]*X[[1]]+v[[2]]*X[[2]],X,False],Norm[#]<r&];
tpy=Select[tangentrealpoints[f,v[[2]]*X[[1]]-v[[1]]*X[[2]],X,False],Norm[#]<r&];
ip=Select[realinflectionpoints[f,X],Norm[#]<r&];
cp=DeleteDuplicates[Join[tpx,tpy,ip]];
cp=Select[cp,   !almostMemberQ[sp,#,.022]&];
(*  Handle singular points *)
nsp=Length[sp];
If[nsp>0,
Do[q=sp[[i]]; 
sol2=Cases[X/.NSolve[{f,(X[[1]]-q[[1]])^2+(X[[2]]-q[[2]])^2-del^2}],{_Real,_Real}];
If[Length[sol2]==0,Print["Potential isolated point at ",q],
Do[stv=GFrame[f,w,X][[1]]; 
If[Norm[w+del*stv-q]>del,gbp=Append[gbp,w], gep=Append[gep,w]],{w,sol2}]
],{i,nsp}]];
 (*  end singular points *)
{gbp,cp,gep,Chop[sp]}]

nline[f_,p_,X_]:=Module[{nl},
nl=tline[f,p,X]/.Thread[X->{X[[2]],-X[[1]]}];
nl-(nl/.Thread[X->p] )]

Curvetrace[ff_,p_,tv_,del_,X_, n_]:= Module[{f,p1,p2,w,df,np,tp,nl,L},
Clear[u137,v137];
f=ff/.Thread[X->{u137,v137}];
p2=p;
L={p2};
Do[p1=p2;
df={D[f,u137],D[f,v137]}/.Thread[{u137,v137}->p1];
np=Expand[df/Norm[df]];
tp={np[[2]],-np[[1]]};
If[tp.tv<0, tp=-tp];
w=p1+del*tp;
nl=tp.{u137,v137};
nl=nl-(nl/.Thread[{u137,v137}->w]);
p2={u137,v137}/.FindRoot[{f,nl}, {{u137,w[[1]]},{v137,w[[2]]}}];
L=Append[L,p2],
{k,n-1}];
L
]

gtv[f_,pp_,AP_,X_]:=Module[{p,q,tv,nv,,ta,tl,tlq,del,vap,tq,nq,k,VAP,nl,j,ct,L,Ans,Q},
del= .2;
p=pp;
Q=False;
{tv,nv}=GFrame[f,p,X];
nl=nline[f,p,X];
tl=tline[f,p,X];
If[(nl/.Thread[X->p+del*tv])<0,nl=-nl];
VAP=Select[AP, ((nl/.Thread[X->#])>0.001) &];   

vap=Length[VAP];
If[vap==0,Print["error 1 at ",p];Return[{p,p}]];
If[vap==1, Return[{p,VAP[[1]]}]];
VAP=SortBy[VAP,Norm[p-#]&];
If[Q,Print["First VAP",VAP]];
k=Min[Floor[Norm[VAP[[1]]-p]/del],10]+1;
If[k==1,If[Q,Print["close enough"]];Return[{p,VAP[[1]]}]];
If[Q,Print["k =",k]];
q=p;
L={p};
j=0;
While[vap>1 && j<10,
ct=Curvetrace[f,q,tv, del, X, k];
If[Norm[ct[[k]]-q]>(k+.5)*del, Print["error 3 at ",p, " -> ",q];Return[L]];
q=ct[[k]];
{tq,nq}=GFrame[f,q,X];
If[tv.tq<-del/2, Print["error 2 at",q]; Return[L]];
L=Append[L,q];
If[Q,Print["tl(q)",tl/.Thread[X->q]]];
If[(tl/.Thread[X->q])<0,tl=-tl];
nl=nline[f,q,X];
If[(nl/.Thread[X->q+del*tq])<0,nl=-nl];
VAP=Select[VAP,((nl/.Thread[X->#])>0)&];

If[Q,Print["VAP2 unsorted ",VAP]];
vap=Length[VAP];
If[vap==0,Print["error 1a ",p," -> ",q]; Return[L]];
If[vap==1,L=Append[L,VAP[[1]]]; Return[L]];
VAP=SortBy[VAP,Norm[q-#]&];
k=Min[Floor[Norm[VAP[[1]]-q]/del],10]+1;
If[Q,Print["k ",k]];
If[k==1,If[Q,Print["close enough"]];L=Append[L,VAP[[1]]];Return[L]];
j++];
If[   !almostMemberQ[AP,L[[-1]],del],Print["Possible error at ",p]];
L
]

QData[f_,r_,X_,opt_:0]:= Module[{L,AP,BP,CP,EP,IP,SP,v,plt0,plt1,plt2,q,nsp, sol,del,stv,n},
n=RandomInteger[{101,999}];
If[Length[opt]==2&&Depth[opt]==2,gvv[n]=Normalize[N[opt]],gvv[n]=Normalize[RandomReal[{-1,1},2]]];
Print["gvv[",n,"] = ",gvv[n]];
plt0:=ContourPlot[x^2+y^2==r^2,{x,-r,r},{y,-r,r},ContourStyle->Red];
{BP,CP,EP,SP}=GCP[f,r,gvv[n],X];
IP=Join[BP,CP];
AP=Join[CP,EP];
L={};
Do[L=Append[L, gtv[f,IP[[i]],AP,X]],{i,Length[IP]}];

nsp=Length[SP];
del= .2; 
If[nsp>0,
Do[q=SP[[i]]; 
sol=Cases[X/.NSolve[{f,(X[[1]]-q[[1]])^2+(X[[2]]-q[[2]])^2-del^2}],{_Real,_Real}];
If[Length[sol]==0,Continue,
Do[L=Append[L,{w,q}],{w,sol}]],{i,nsp}]];
L];

QPlot[f_,r_,X_,opt_:0]:= GListPlot[f,QData[f,r,X,opt],Black];

gworks[f_,r_,X_,v_]:= Module[{L,AP,BP,CP,EP,IP,SP,plt1,plt2,q,nsp, sol,del,stv},

{BP,CP,EP,SP}=GCP[f,r,v,X];
IP=Join[BP,CP];
AP=Join[CP,EP];
L={};
Do[L=Append[L, gtv[f,IP[[i]],AP,X]],{i,Length[IP]}];

nsp=Length[SP];
del= .2; 
If[nsp>0,
Do[q=SP[[i]]; 
sol=Cases[X/.NSolve[{f,(X[[1]]-q[[1]])^2+(X[[2]]-q[[2]])^2-del^2}],{_Real,_Real}];
If[Length[sol]==0,Continue,
Do[L=Append[L,{w,q}],{w,sol}]],{i,nsp}]];

{L,BP,EP,SP}]

Gsegs[f_,rr_,X_,v_]:=Module[{r,g,Lt,BP,EP,SP,i,t,s,sb,eb,lbp,LL,pos,SL,del,q,w,lip},
r=rr;
{Lt,BP,EP,SP}=gworks[f,r,X,v];

del= .2;
SL={};
Do[
sb=Select[Lt,(Norm[#[[2]]-q]<.001&& almostMemberQ[EP,#[[1]],.001])&][[All,1]];
eb=Select[Lt,(Norm[#[[2]]-q]<.001&& almostMemberQ[BP,#[[1]],.001])&][[All,1]];


g=(X[[1]]-q[[1]])^2+(X[[2]]-q[[2]])^2-.04;

Do[
w= .1*Normalize[Total[GFrame[g,sb[[i]],X]]];
If[(f/.Thread[X->sb[[i]]+w])>0,
s=gtv[g,sb[[i]],eb,X],
s=gtv[-g,sb[[i]],eb,X]];
SL=Append[SL,Insert[s,q,2]],{i,Length[sb]}],
{q,SP}];
If[Length[SP]>0,
Lt=Select[Lt,   !almostMemberQ[SP,#[[2]],.001]&];
Lt=Join[Lt,SL];

BP=Select[BP,   !almostMemberQ[eb,#,.001]&];
EP=Select[EP,   !almostMemberQ[sb,#,.001]&]

];

lbp=Length[BP];
i=1;
While[Length[Lt]>0 && i<11,

LL[i]=Lt[[1]];

Lt=Drop[Lt,1];
s=Last[LL[i]];

If[MemberQ[EP,s],i++;  Continue[]]; 
pos=npos[Flatten[Take[Lt,All,1],1],s,.001];
If[Length[pos]>1,Print["pos error 1"];Abort[]];
While[Length[pos]>0,
 
LL[i]=Join[LL[i],Extract[Lt,pos[[1]]]];

Lt=Drop[Lt,pos[[1]]];

s=Last[LL[i]];

If[MemberQ[EP,s],  Break[]]; 
pos=npos[Flatten[Take[Lt,All,1],1],s,.001];
If[Length[pos]>1, Print["pos error 2"];Abort[]];

]; 
i++];
lip=Length[infiniterealpoints[f,X]];
Print[lip," infinite points ",i-1, " observed components ",Length[SP]," Singular points " ];
If[lip>i-1,Print["Possible Missing Component"]];
Table[LL[j],{j,i-1}]
]

GData[f_,r_,X_,v_]:=Gsegs[f,r,X,Normalize[v]];

GPlot[f_,r_,X_,opt_:0]:=Module[{L,plt1,plt2,plt0,v,n},
n=RandomInteger[{101,999}];
If[Length[opt]==2&&Depth[opt]==2,gvv[n]=Normalize[N[opt]],gvv[n]=Normalize[RandomReal[{-1,1},2]]];
Print["gvv[",n,"] = ",gvv[n]];
L=Gsegs[f,r,X,gvv[n]];
plt0 := ContourPlot[x^2+y^2==r^2,{x,-r,r},{y,-r,r}, ContourStyle->Red];
plt1:=RegionPlot[f>0,{x,-r,r},{y,-r,r}, BoundaryStyle->None];
plt2:=ListLinePlot[L,PlotRange->All, PlotStyle->Thick];
Show[plt0,plt1,plt2]
]

GListPlot[f_,L_,clr_:Black]:=Module[{r,plt0,plt1,plt2},
r=Max[Norm/@ Flatten[L,1]]*1.05;
plt0:=ContourPlot[x^2+y^2==r^2/1.05, {x,-r,r},{y,-r,r}, ContourStyle->Pink];
plt1 := RegionPlot[f>0,{x,-r,r},{y,-r,r},BoundaryStyle->None];
plt2 := ListLinePlot[L, PlotRange->All,PlotStyle->clr];
Show[plt0,plt1,plt2]
];

RefineComponent[f_,L_,X_]:=Module[{s,p,q,LL,np,j},
LL=L;
Do[
p=closestrealpoint[f,(L[[i]]+L[[i+1]])/2,X];
If[Norm[p-L[[i+1]]]<.01,Continue[]];
q=Last[gtv[f,p,{L[[i+1]]},X]];
np=npos[LL,q,.001];
If[Length[np]>0, j=np[[1,1]];
LL=Insert[LL,p,j]],
{i,Length[L]-1}];
LL];

(*  Finds component(s) containing point p and adds p to list if not already there *)
FindComp[f_,p_,L_,X_]:=Module[{r,FL,q,T,S,np,nl,dl,i,j,K},
Check[
nl=Length[L];
dl=Depth[L];
FL=If[dl==4,Flatten[L,1],L];
q=Last[gtv[f,p,FL,X]];
i=0;
While[i≤ nl,
np=npos[L[[i]],q,.001];
If[Length[np]>0, j=np[[1,1]];
Return[{i,Insert[L[[i]],p,j]}]];
i++],
Print["error"]]
];