- /
-
Orbiting Roots
on 13 Nov 2023
- 23
- 176
- 2
- 5
- 1067
drawframe(1);
function drawframe(f)
t = interp1([0 48],[0 2*pi],f);
p1 = [1 0 0 0 0 1 0];
p2 = [1 -0.5*(cos(t) + 1i*sin(t))];
p = conv(p1,p2);
lims = 3*[-1 1];
n = 300;
method = "halley";
[zrt,it,xs,ys] = root_basins(p,lims,lims,n,method);
imagesc(it)
set(gca,XTick=[],YTick=[])
colormap(jet)
set(gca,CLim=[1 10])
axis square
end
function [rootMap,iterMap,xs,ys] = root_basins(p,xlims,ylims,n,method)
rt = roots(p);
% Take the necessary derivatives
dp = polyder(p);
ddp = polyder(dp);
% Make some helpful local functions
f = @(q) polyval(p,q);
df = @(q) polyval(dp,q);
ddf = @(q) polyval(ddp,q);
xs = linspace(xlims(1),xlims(2),n);
ys = linspace(ylims(1),ylims(2),n);
[X,Y] = meshgrid(xs,ys);
z = X + 1i*Y;
tol = 1e-6;
iterMap = zeros(size(z)); % Track iteration count
rootMap = zeros(size(z)); % Track roots
freezeMap = false(size(z)); % Track where you're done
for k = 1:100
% Only calculate on unfrozen regions
zz = z(~freezeMap);
if method == "newton"
zz = zz - f(zz)./df(zz);
elseif method == "halley"
zz = zz - 2*f(zz).*df(zz)./(2*(df(zz).^2 - f(zz).*ddf(zz)));
else
disp("No such method")
end
% Put the newly calculated values back into z
z(~freezeMap) = zz;
% Loop over the roots to see where we're close enough to be done
for m = 1:length(rt)
zerosFoundMap = abs(rt(m)-z) < tol;
% Only update when you find zeros in unfrozen regions
zerosFoundMap = zerosFoundMap & ~freezeMap;
if any(zerosFoundMap(:))
rootMap(zerosFoundMap) = m;
iterMap(zerosFoundMap) = k;
freezeMap = freezeMap | zerosFoundMap;
end
end
end
end
