Date:         Fri, 12 Jul 2002 14:14:48 -0700
Reply-To:     Clayton Wells <clay_wells@HOTMAIL.COM>
Sender:       "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From:         Clayton Wells <clay_wells@HOTMAIL.COM>
Organization: http://groups.google.com/
Subject:      SAS IML prog to generate Newton's Method fractal
Content-Type: text/plain; charset=ISO-8859-1

I came across this while cleaning out some old stuff.  It's something
I wrote several years ago.  Thought I would share it.  Kinda fun to
play around with - although to get a good image, a fast machine helps.
 IML generates the points and GPLOT is used for the plot of course.

This is Newton's method applied to the function f(x)=x^3-1.

All you must change is the libname and right below that the variable
caled "ps".  this is the number of points the program generates.  it's
set at 1000 now - which shows nothing really.  but you can start from
there and work up to 1,000,000 or whatever.

enjoy - or not.

cw

 ==========================================================================

libname lib 'c:\';     /* assign an appropriate path for the libname
here  */


proc iml symsize=3000;

ps=1000;  /* number of points to be generated  */



/* NOTHING BELOW THIS LINE NEEDS MODIFYING */


a=repeat(0,15,2);
b=99;
c=99;
d=repeat(0,ps,1);
e=repeat(0,ps,1); /*  start points  */
f=repeat(0,ps,1);
gg=repeat(0,ps,1); /* root */
hh=repeat(0,ps,1); /* root */
g=repeat(0,15,2);
gp=repeat(0,15,2);


do j=1 to ps;

size1=3;
size2=3;
do until (abs(size1)<1 & abs(size2)<1);        /* AREA CONTROLLED HERE
 */
    a[1,1]=rannor(0);
    a[1,2]=rannor(0);
    size1=a[1,1];
    size2=a[1,2];
end;

b=99;
c=99;
e[j]=a[1,1];
f[j]=a[1,2];
g[1,1]=a[1,1]**3-3*a[1,1]*a[1,2]**2-1;
g[1,2]=3*a[1,1]**2*a[1,2]-a[1,2]**3;
gp[1,1]=3*(a[1,1]**2-a[1,2]**2);
gp[1,2]=6*a[1,1]*a[1,2];

i=2;
do while (abs(b)>0.0001 & abs(c)>0.0001);

  a[i,1]=a[i-1,1] - ((a[i-1,1]**3-3*a[i-1,1]*a[i-1,2]**2-1)*(3*(a[i-1,1]**2-a[i-1,2]**2))
+
          (3*a[i-1,1]**2*a[i-1,2]-a[i-1,2]**3)*(6*a[i-1,1]*a[i-1,2]))/
          (((3*a[i-1,1]**2-3*a[i-1,2]**2)**2)+(6*a[i-1,1]*a[i-1,2])**2);

  a[i,2]=a[i-1,2] - ((3*a[i-1,1]**2*a[i-1,2]-a[i-1,2]**3)*(3*(a[i-1,1]**2-a[i-1,2]**2))
-
          (a[i-1,1]**3-3*a[i-1,1]*a[i-1,2]**2-1)*(6*a[i-1,1]*a[i-1,2]))/
          (((3*a[i-1,1]**2-3*a[i-1,2]**2)**2)+(6*a[i-1,1]*a[i-1,2])**2);

  g[i,1]=a[i-1,1]**3-3*a[i-1,1]*a[i-1,2]**2-1;
  g[i,2]=3*a[i-1,1]**2*a[i-1,2]-a[i-1,2]**3;
  gp[i,1]=3*(a[i-1,1]**2-a[i-1,2]**2);
  gp[i,2]=6*a[i-1,1]*a[i-1,2];

b=g[i,1];  c=g[i,2];
d[j]=i-1;
if d[j]=14 then b=.00000000001;
i=i+1;

gg[j]=a[d[j]+1,1];
hh[j]=a[d[j]+1,2];
end;


end;




dset=repeat(0,ps,5);
dset=e || f || gg || hh || d;


create t1 from dset[colname='start_re'];
append from dset;



* print dset;
quit;

data lib.new; set t1;
run;


/*   ******   */

data temp; set lib.new;
run;

data temp2(rename=(col2=start_im col3=root_re col4=root_im
col5=iterate)); set temp;
run;


data temp2; set temp2;
if (0.9 <= root_re <= 1.01 & -0.01 <= root_im <= 0.01) then conv=1;
if (-0.51 <= root_re <=-0.49 & 0.85 <= root_im <=0.87) then conv=2;
if (-0.51 <= root_re <=-0.49 & -0.87 <= root_im <= -0.85) then conv=3;
run;


proc freq data=temp2;
tables conv;
run;

/* ************** */


data a; set temp2;
if conv=1 then do;
    st_re_1=start_re;
    st_im_1=start_im;
end;
if conv=2 then do;
    st_re_2=start_re;
    st_im_2=start_im;
end;
if conv=3 then do;
    st_re_3=start_re;
    st_im_3=start_im;
end;
run;


data lib.a; set a;
run;



symbol1 color=CYAN value=dot height=.01;
symbol2 color=yellow value=dot height=.01;
symbol3 color=PINK value=dot height=.01;

axis1 order=(-1 to 1)
  label=none
  MINOR=NONE;
axis2 order=(-1 to 1)
  label=none
  MINOR=NONE;


proc gplot data=a;
  plot st_im_1*st_re_1 st_im_2*st_re_2 st_im_3*st_re_3
     / overlay haxis=axis1 vaxis=axis2 cframe=purple
       HREF=0 VREF=0;
run;

quit;





/* *************** */