% scripts for October 12, Newton's method, Julia sets, Mandelbrot set %--------------------------------------------------------- % function to apply Newton's method to a polynomial % derivative is computed by this function function [y,k] = poly_newton(p,x) p_prime=polyder(p); k = 0; xprev=x+1; while abs(x - xprev) > eps*abs(x) xprev = x; x = x - polyval(p,x)/polyval(p_prime,x); k = k + 1; if k>50 % stop if doesn't converge y=NaN; return; end % if end % while y=x; return; % try p=[1,0,0,1] cube roots of 1 % try p=[2,8,2,0,3,1] %-------------------------------------------------------- % function to plot result of Newton's method delta=abs(2*a)/n; % save_root=zeros(n); for k=1:n for j=1:n z=(-abs(a)+2*abs(a)*j/n)+1i*(-abs(a)+2*abs(a)*k/n); [y,iter]=poly_newton(p,z); if isnan(y) save_root(k,j)=0; else [m,root_num]=min(abs(all_roots-y)); save_root(k,j)=root_num; end end k end all_roots; save_root; figure; title('Domains of attraction for Newtons method', p) image(save_root,'CDataMApping','scaled') %axis image %colormap(flipud(jet(length(all_roots)))) toc; return; %-------------------------------------- % script to test preallocation tic clear t N=10000000; % t=zeros(N,1); t(1)=0; counter=1; while counter < N t(counter+1)=t(counter)+1; counter=counter+1; end toc %--------------------------------------- %function to draw Mandelbrot set with for loops function mandel1(center, radius,n) tic; cx=real(center); cy=imag(center); delta=2*radius/n; x=[cx-radius:delta:cx+radius]; y=[cy-radius:delta:cy+radius]; Lx=length(x); Ly=length(y); %c=zeros(Ly,Lx); depth=50; for j=1:Lx for m=1:Ly z0=x(j)+i*y(m); z=0; c(m,j)=depth; for k=1:depth z=z^2+z0; if abs(z)>2 c(m,j)=k; break; end % if end % for k end % for m end % for j image(c) axis image colormap(flipud(jet(depth))) toc; return; %----------------------------------------------- %-------------------------------- %function to draw Mandelbrot set using vectorization function mandel2(center, radius,n) tic; cx=real(center); cy=imag(center); delta=2*radius/n; x=[cx-radius:delta:cx+radius]; y=[cy-radius:delta:cy+radius]; [X,Y] = meshgrid(x,y); z0=X+i*Y; z=zeros(n+1,n+1); c=zeros(n+1,n+1); depth=50; for k=1:depth z=z.^2+z0; c(abs(z)<2)=k; end image(c) axis image colormap(flipud(jet(depth))) toc; return; %--------------------------------------------- % function to draw Mandelbrot set with GUI controls function mandel(varargin) if nargin == 0 center=-.5; radius=2; n=500; elseif nargin==2 center=varargin{1}; radius=varargin{2}; n=500; elseif nargin==3 center=varargin{1}; radius=varargin{2}; n=varargin{3}; else return; end tic; depth=50; cx=real(center); cy=imag(center); delta=2*radius/n; x1=cx-radius; x2=cx+radius; y1=cy-radius; y2=cy+radius; flag =1; while flag==1 delta=(x2-x1)/n x=[x1:delta:x2]; y=[y2:-delta:y1]; Lx=length(x); Ly=length(y); [X,Y] = meshgrid(x,y); z0=X+i*Y; size(z0); z=zeros(Ly,Lx); c=zeros(Ly,Lx); for k=1:depth z=z.^2+z0; c(abs(z)<2)=k; end m=min(min(c)); M=max(max(c)); d=(M/(M-m))*(c-m); image(d) % choose how to label axes axistype=1; switch axistype case 0 % no labels axis off case 1 % label with planar coordinates horz=x1:(x2-x1)/4:x2; for k=1:length(horz) xlabels{k}=num2str(horz(k)); end Lx; xt=[[1:Lx/4:Lx],Lx]; xticks(xt) xticklabels(xlabels) vert=y2:(y1-y2)/4:y1; for k=1:length(vert) ylabels{k}=num2str(vert(k)); end Ly; yt=[[1:Ly/4:Ly],Ly]; yticks(yt) yticklabels(ylabels) case 2 % label with matrix coordiates axis image end % switch % set colors to use setc=flipud(jet(depth)); %colormap(jet(depth)) %colormap(hsv(depth)) %colormap(gray(depth)) %colormap(flipud(winter(depth))) s=size(setc,1); setc([s-1,s],:)=zeros(2,3); colormap(setc) % add a legend saying how colors give iterates colorbar toc; % get data from mouse to draw new picture [a(1),b(1),button]=ginput(1) % if middle mouse button is hit, input new depth depth=max(depth,min(min(c))+100); switch button case 1 [a(2),b(2)]=ginput(1); case 2 depth=input('input new depth = '); n=input('input new pixels = '); case 3 z1=x1+delta*a(1); z2=y2-delta*b(1); fig=gcf; julia(z1,z2); figure(fig) [a(1),b(1),button]=ginput(1) [a(2),b(2)]=ginput(1); end % switch if button ~= 2 %[a,b]=ginput(2); a=sort(a); b=sort(b); a2=x1+delta*a(2); a1=x1+delta*a(1); b1=y2-delta*b(2); b2=y2-delta*b(1); cx=(a1+a2)/2; cy=(b1+b2)/2; r=abs(a2-a1)/2; x1=cx-r; x2=cx+r; y1=cy-r; y2=cy+r; end % if button end % while return; %------------------------------------ %function to draw quadratic Julia sets function julia(a,b) tic; depth=100; x1=-2; x2=2; y1=-2; y2=2; delta=4/1000; x=x1:delta:x2; y=y2:-delta:y1; Lx=length(x); Ly=length(y); [X,Y] = meshgrid(x,y); z0=X+i*Y; size(z0); z=z0; c=zeros(Ly,Lx); for k=1:depth z=z.^2+(a+i*b); c(abs(z)<2)=k; end m=min(min(c)) M=max(max(c)) d=c; %(M/(M-m))*(c-m); figure; title('Julia set'); image(d) % choose how to label axes axistype=1; switch axistype case 0 % no labels axis off case 1 % label with planar coordinates horz=x1:(x2-x1)/4:x2; for k=1:length(horz) xlabels{k}=num2str(horz(k)); end Lx; xt=[[1:Lx/4:Lx],Lx]; xticks(xt) xticklabels(xlabels) vert=y2:(y1-y2)/4:y1; for k=1:length(vert) ylabels{k}=num2str(vert(k)); end Ly; yt=[[1:Ly/4:Ly],Ly]; yticks(yt) yticklabels(ylabels) case 2 % label with matrix coordiates axis image end % switch % set colors to use setc=flipud(jet(depth)); %colormap(jet(depth)) %colormap(hsv(depth)) %colormap(gray(depth)) %colormap(flipud(winter(depth))) s=size(setc,1) setc(s,:)=zeros(1,3); colormap(setc) % add a legend saying how colors give iterates colorbar toc; return