Date: Fri, 12 Jul 2002 14:14:48 -0700 Reply-To: Clayton Wells Sender: "SAS(r) Discussion" From: Clayton Wells 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; /* *************** */

Back to: Top of message | Previous page | Main SAS-L page