% 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