Durand–Kerner method matlab code help? Anyone can help? Can't figure out what is wrong.
4 views (last 30 days)
Show older comments
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;
5 Comments
Answers (1)
Andrei Bobrov
on 10 Oct 2011
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
0 Comments
See Also
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!