% This program demonstrates how all complex roots of polynomial p1 % can be found by homotopy from polynomial p0 with known roots. % This program trying to predict what happens when two roots z1 and z2 collide, % and jump over the dangereous value T of t, just rotating two close roots % around their midpoint zmid by 90 degrees function value=main() % p0 is a polynomial with known roots; coefficients are stored in ascending order, % the first one is the constant term, the last one correspond to the maximal degree p0 = [0 1 0 1 0.2]; n = length(p0); deg=n-1; p=p0; % an intermediate polynomial of the same degree as p0 roots=[0 i -i -5]; % roots of p0 newroots=roots; Roots=roots; % Here are our initial data: polynomial p1 and number of steps p1 = [2 -3 -2 1 1]; % for this initial data there are two collisions number_of_steps = 100; dt=1./number_of_steps; t=0; while t<1-0.5*dt, t=t+dt; for x = 1:n, p(x) = (1-t)*p0(x) + t*p1(x); end for k=1:deg, z=roots(k); dz=P(p,z)/DP(p,z); newroots(k)=z-dz; end % ------- let's check if there is a collision and which roots are colliding distmin= 10^15; kmin=0; mmin=0; dzmin=0; for k=1:deg-1, for m=k+1:deg, dz=newroots(k)-newroots(m); dist=abs(dz); if dist=0.75*abs(olddz), roots=newroots; Roots=[Roots; roots]; end if abs(newdz)< 0.75*abs(olddz), t=t-dt; for x = 1:n,p(x) = (1-t)*p0(x) + t*p1(x); end % one step back z1=roots(kmin) z2=roots(mmin) Z1=newroots(kmin) Z2=newroots(mmin) pause for l=1:2, z1=z1-P(p,z1)/DP(p,z1); z2=z2-P(p,z2)/DP(p,z2); end d1=(z1-z2)^2 t1=t t=t+dt/2; for x = 1:n, p(x)=(1-t)*p0(x)+t*p1(x); end for l=1:2, z1=z1-P(p,z1)/DP(p,z1); z2=z2-P(p,z2)/DP(p,z2); end d2=(z1-z2)^2 t2=t % near to the moment of collision T, expression (z1-z2)^2 is % a linear function in T-t, let T = t2+x, then x=(t2-t1)*d2/(d1-d2); x=abs(x); T=t2+x t3=t2+2*x zmid=(z1+z2)/2 dz=(z1-z2)/2 newz1=zmid+i*dz newz2=zmid-i*dz pause t=t3; dt=1./(1-t3)/number_of_steps; for x = 1:n, p(x) = (1-t)*p0(x) + t*p1(x); end roots(kmin)=newz1; roots(mmin)=newz2; for l=1:2, for k=1:deg,z=roots(k);roots(k)=z-P(p,z)/DP(p,z);end Roots=[Roots; roots]; end end%if end%while % final improvement of roots using 5 steps of Newton for l=1:5, for k=1:deg, z=roots(k); roots(k)=z-P(p1,z)/DP(p1,z); end Roots=[Roots; roots]; end roots=roots % output %-------------------------------------------------------- [d,dd]=size(Roots);x=1:d;y=1:d; for k=1:deg, for l=1:d, x(l)=real( Roots(l,k) ); y(l)=imag( Roots(l,k) ); end if k==2 hold on end plot(x,y,'o-'); end r=2; axis( [-r,r,-r,r] ); grid; hold off; %-------------------------------------------------------- function value = P(a,z) n= length( a ); value = a(n); for i = n-1:-1:1, value = value*z; value = value + a(i); end function value = DP(a,z) n = length( a ); value = (n-1)*a(n); for i = n-1:-1:2, value = value*z; value = value + (i-1)*a(i); end