Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Aug 2016
Lecture topic: Roots of Unity and Roots of Complex Numbers

Newton's Method & Roots of Unity: (z^n)-1=0

The n roots of a complex number y  are given by the following, where 0≤k≤n-1.

y:=c0+c1*I;
r:=Sqrt[(c0^2)+(c1^2)];
θ:=ArcTan[(c1/c0)];

root[k_]:=(r^(1/n))*(Cos[((θ+2Pi*k)/n)]+I*Sin[((θ+2Pi*k)/n)])

Use Newton's Method to determine which root the intial points go to, and partition the phase space on this criteria (n=5)

The number of roots is given by n.

Clear[c0,c1,n]
npts=20000;
c0=1.;
c1=0.;
n=5.;
f[x_]:=(x^n)-y;
fp[x_]:=f'[x];
Do[set[i]={},{i,0,n-1}]
acc=4;
tol=(10.^(-acc));
max=500;
Do[Clear[p0,p1,i,p,q,a,b];
a=RandomReal[{-1,1}];
b=RandomReal[{-1,1}];
p0=a+b I;i=1;
While[i≤max,
p=p0-((f[p0])/(fp[p0]));If[Abs[f[p]]<tol,{Do[If[Abs[p-root[k]]<(10^(-acc+1)),set[k]=Append[set[k],{a,b}]],{k,0,n-1}],Break[]}];p0=p;i=i+1],{jj,1,npts}]

Plot the phase space, colouring each point depending on which root it goes to  (n=5)

setlist={};
Do[setlist=Append[setlist,set[i]],{i,0,n-1}]
huelist={};
Do[huelist=Append[huelist,Hue[(i/n)]],{i,0,n-1}]
ListPlot[setlist,AspectRatio->Automatic,Axes->True,Frame->True,PlotStyle->huelist,Background->GrayLevel[0],ImageSize->{500,500}]

Here is a plot with a n=4

Clear[c0,c1,n]
npts=20000;
c0=1.;
c1=0.;
n=4.;
f[x_]:=(x^n)-y;
fp[x_]:=f'[x];
Do[set[i]={},{i,0,n-1}]
acc=3;
tol=(10.^(-acc));
max=500;
Do[Clear[p0,p1,i,p,q,a,b];
a=RandomReal[{-1,1}];
b=RandomReal[{-1,1}];
p0=a+b I;
i=1;
While[i≤max,
p=p0-((f[p0])/(fp[p0]));If[Abs[f[p]]<tol,{Do[If[Abs[p-root[k]]<(10^(-acc+1)),set[k]=Append[set[k],{a,b}]],{k,0,n-1}],Break[]}];p0=p;i=i+1],{jj,1,npts}]

setlist={};
Do[setlist=Append[setlist,set[i]],{i,0,n-1}]
huelist={};
Do[huelist=Append[huelist,Hue[(i/n)]],{i,0,n-1}]
ListPlot[setlist,AspectRatio->Automatic,Axes->True,Frame->True,PlotStyle->huelist,Background->GrayLevel[0],ImageSize->{500,500}]

Newton's Method & Roots of a Complex Number: (z^n)-y=0

Clear[c0,c1,n]
npts=20000;
c0=(1/2.);
c1=(1/3.);
n=3.;
f[x_]:=(x^n)-y;
fp[x_]:=f'[x];
Do[set[i]={},{i,0,n-1}]
acc=4;
tol=(10.^(-acc));
max=500;
Do[Clear[p0,p1,i,p,q,a,b];
a=RandomReal[{-1,1}];
b=RandomReal[{-1,1}];
p0=a+b I;
i=1;
While[i≤max,
p=p0-((f[p0])/(fp[p0]));If[Abs[f[p]]<tol,{Do[If[Abs[p-root[k]]<(10^(-acc+1)),set[k]=Append[set[k],{a,b}]],{k,0,n-1}],Break[]}];p0=p;
i=i+1],{jj,1,npts}]

setlist={};
Do[setlist=Append[setlist,set[i]],{i,0,n-1}]
huelist={};
Do[huelist=Append[huelist,Hue[(i/n)]],{i,0,n-1}]
ListPlot[setlist,AspectRatio->Automatic,Axes->True,Frame->True,PlotStyle->huelist,Background->GrayLevel[0],ImageSize->{500,500}]

A Mathematica Module

plotrooty[c0_,c1_,n_,npts_]:=Module[{y,r,f,fp},
y=c0+c1 I;
r=(Sqrt[(c0^2)+(c1^2)]);
theta=ArcTan[(c1/c0)];
root[i_]:=(r^(1/n)) (Cos[((theta+2 π i)/n)]+I Sin[((theta+2 π i)/n)]);
f[x_]:=(x^n)-y;
fp[x_]:=f'[x];
Do[set[i]={},{i,0,n-1}];
acc=4;
tol=(10.^(-acc));
max=500;
Do[Clear[p0,p1,i,p,q,a,b];
a=RandomReal[{-1,1}];
b=RandomReal[{-1,1}];
p0=a+b I;
i=1;
While[i≤max,p=p0-((f[p0])/(fp[p0]));If[Abs[f[p]]<tol,{Do[If[Abs[p-root[k]]<(10^(-acc+1)),set[k]=Append[set[k],{a,b}]],{k,0,n-1}],Break[]}];p0=p;
i=i+1],{jj,1,npts}];
setlist={};
Do[setlist=Append[setlist,set[i]],{i,0,n-1}] ;
huelist={};
Do[huelist=Append[huelist,Hue[(i/n)]],{i,0,n-1}];ListPlot[setlist,AspectRatio->1,Axes->True,Frame->True,PlotRange->{{-1,1},{-1,1}},PlotStyle->huelist,Background->GrayLevel[0],ImageSize->{500,500}]]

plotrooty[1.0,0.0,4.0,10000]

A Mathematica Module for First Quadrant Only

plotrooty2[c0_,c1_,n_,npts_]:=Module[{y,r,f,fp},
y=c0+c1 I;
r=(Sqrt[(c0^2)+(c1^2)]);
theta=ArcTan[(c1/c0)];
root[i_]:=(r^(1/n)) (Cos[((theta+2 π i)/n)]+I Sin[((theta+2 π i)/n)]);
f[x_]:=(x^n)-y;
fp[x_]:=f'[x];
Do[set[i]={},{i,0,n-1}];
acc=4;
tol=(10.^(-acc));
max=500;
Do[Clear[p0,p1,i,p,q,a,b];
a=RandomReal[{0,1}];
b=RandomReal[{0,1}];
p0=a+b I;
i=1;
While[i≤max,
p=p0-((f[p0])/(fp[p0]));If[Abs[f[p]]<tol,{Do[If[Abs[p-root[k]]<(10^(-acc+1)),set[k]=Append[set[k],{a,b}]],{k,0,n-1}],Break[]}];p0=p;
i=i+1],
{jj,1,npts}];
setlist={};
Do[setlist=Append[setlist,set[i]],{i,0,n-1}] ;
huelist={};
Do[huelist=Append[huelist,Hue[(i/n)]],{i,0,n-1}];ListPlot[setlist,AspectRatio->1,Axes->True,Frame->True,PlotRange->{{0,1},{0,1}},PlotStyle->huelist,Background->GrayLevel[0],ImageSize->{500,500}]]

plot1=plotrooty2[1.0,0.0,4.0,20000]

A Mathematica Module using Halton Sequence

This module just works in the first quadrant, and uses the Halton quasi-random sequence. We will look at Corput and Halton sequence when we do Monte Carlo integration, but it produces such nice pictures here that it's worth including it now. Basic Idea is--rather than true random numbers, use a quasi-random sequence that samples space in a regular manner.

plotrootyhalton[c0_,c1_,n_,npts_]:=Module[{y,r,f,fp,Corput},
y=c0+c1 I;
r=(Sqrt[(c0^2)+(c1^2)]);
theta=ArcTan[(c1/c0)];
root[i_]:=(r^(1/n)) (Cos[((theta+2 π i)/n)]+I Sin[((theta+2 π i)/n)]);Corput[nn_,b_]:=N[IntegerDigits[nn,b].(b^(Range[-Floor[Log[b,nn]+1],-1]))];
f[x_]:=(x^n)-y;
fp[x_]:=f'[x];
Do[set[i]={},{i,0,n-1}];
acc=4;
tol=(10.^(-acc));
max=500;
Do[Clear[p0,p1,i,p,q,a,b];
a=Corput[jj,2];
b=Corput[jj,3];
p0=a+b I;
i=1;
While[i≤max,
p=p0-((f[p0])/(fp[p0]));If[Abs[f[p]]<tol,{Do[If[Abs[p-root[k]]<(10^(-acc+1)),set[k]=Append[set[k],{a,b}]],{k,0,n-1}],Break[]}];p0=p;
i=i+1],
{jj,1,npts}];
setlist={};
Do[setlist=Append[setlist,set[i]],{i,0,n-1}] ;
huelist={};
Do[huelist=Append[huelist,Hue[(i/n)]],{i,0,n-1}];ListPlot[setlist,AspectRatio->1,Axes->True,Frame->True,PlotRange->{{0,1},{0,1}},PlotStyle->huelist,Background->GrayLevel[0],ImageSize->{500,500}]]

plot2=plotrootyhalton[1.0,0.0,4.0,20000]

Comparison of Random vs Halton (quasi-random)

Show[GraphicsRow[{plot1,plot2}],ImageSize->{900,Automatic}]

An Extremely Nice Picture (takes a looooong while! When it's done, drag it to make it bigger.)

plot2=plotrootyhalton[1.0,0.0,4.0,100000]

Halton Sequence to Generate quasi-random numbers (We'll look at this during Monte Carlo Integration)

Random Numbers

SeedRandom[2];
random=Table[RandomReal[{0,1},2],{i,1,2000}];

rplot=ListPlot[random,Frame->True,AspectRatio->1,PlotStyle->Blue,ImageSize->{500,Automatic}]

Quasi-Random: Corput Sequence (one Dimensional)

Corput[n_,b_]:=IntegerDigits[n,b].(b^(Range[-Floor[Log[b,n]+1],-1]))

corput=Table[Corput[i,2],{i,1,20}]

corput=Table[Corput[i,3],{i,1,20}]

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

corput=Table[Corput[i,2],{i,1,20}]

corput=Table[Corput[i,3],{i,1,20}]

Quasi-Random: Halton Sequence (two dimensional)

Halton[n_]:={Corput[n,2],Corput[n,3]}

halton=Table[Halton[k],{k,1,10}]

halton=Table[Halton[k],{k,1,2000}];

hplot=ListPlot[halton,Frame->True,AspectRatio->1,PlotStyle->Red,ImageSize->{500,Automatic}]

Comparison

Show[rplot,hplot]