Durand–Kerner method matlab code help? Anyone can help? Can't figure out what is wrong.
4 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Andrei Bobrov
2011-10-10
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!