Durand–Kerner method matlab code help? Anyone can help? Can't figure out what is wrong.

Hello Experts,

A is a vector of polynomial coefficients and I need to find roots using the Durand-Kerner method. Please have a look at the code, something very strange occurs, I do find roots but not all of them are right one. I have used this: http://ezekiel.vancouver.wsu.edu/~cs330/projects/dk/DK-method.pdf

I have written the following function, please just have a look what is wrong, I am following the algorithm but maybe I miss something.

Be so kind and help me, thanks in advance!

Take for example: Poly coeffs A=[1,-6,11,-6] roots are: 1; 2; 3 After running this function I get: -13.8879 +10.3923i, -13.8879 -10.3923i, 4.1121 - 0.0000i

Already tried to figure out what is wrong but can't understand...

function rvec = roots(A)

rvec = zeros(1);
x=zeros(1);
d = size(A,2)-1;
A = A/A(1); 
Nmax = 50;
tol = 10^(-7);
max_coeff = max(A);
Radius = 1+max_coeff;
theta = 0;

for v=1:d

      theta = (2*pi*v)/d;
      x(v) = Radius*complex(cos(theta),sin(theta));

end

m=ones(1,d);

for t=1:d
    for u=1:d
        if t~=u
            m(t) = m(t)*(x(t)-x(u));
        end
    end
end

stop = 0;

for k=1:Nmax
    if stop == 0;
        delta_x_max = 0;
        for z=1:d
            delta_x = -polyval(A,x(z))/m(z);

              if norm(delta_x,2)>delta_x_max
                  delta_x_max = norm(delta_x,2);
              end

          end
          for t=1:d
              x(t) = x(t) + delta_x;
          end

          if delta_x_max<=tol
              stop=1;
          end
      else
          break;
      end
  end
  rvec = x;

[Walter Roberson]

No errors shout out at me. I do see several opportunities for optimization or code cleanup though. For example, eliminate the "if stop == 0" and at the place you set stop=1, instead have a direct "break" statement.

[Steve]

Hello Walter,

If you will try my algorithm with A = [1,-3,3,-5] you will get real root that is right but the complex one isn't right.

And what other optimization should I make to my code?

[Walter Roberson]

I can't test it this weekend as remote access is still down to my server.

The "for v" loop can be completely vectorized.

theta = 2*pi/d .* (1:d);

x = Radius*complex(cos(theta),sin(theta));

Consider

temp = bsxfun(@minus,x(:),x);

temp = temp + eye(d);

m = prod(temp,2);

[Steve]

Dear Walter, please try to run this function somewhere I tried to understand what is going wrong with this but I get really strange roots.

I have sent to you the full email with code and example.

Please be so kind, have a look.

[Walter Roberson]

I do not have access to a MATLAB instance tonight or tomorrow.

Answer 1

[Andrei Bobrov]

variant

c1 - coefficient of input polynom (eg: f(z) = z^3 - 3*z^2 + 3*z - 5; c1 = [1 -3 3 -5])

c1 = c1(:)/c1(1);
c = c1(end:-1:2);
nmax = 500;
n = numel(c);
z = (1+max(abs(c)))*arrayfun(@(i1)complex(cos(i1),sin(i1)),2*pi/n*(0:n-1)');
tls = 1e-6;
e1 = eye(n);
for j1 = 1:nmax
  dz = polyval(c1, z)./prod(bsxfun(@minus,z,z.')+e1,2);
  z = z - dz;
  if norm(dz,inf) <= tls 
      break
  end
end