求聚根的Laguerre方法(Matlab)
我必须使用Laguerre的方法编写一段代码，以找到poly的真实和复杂的根：
P=X^5-5*X^4-6*X^3+6*X^2-3*X+1
我一点也不怀疑。我在matlab中做了算法，但是5个根中有3个是相同的，我不认为这是正确的。
syms X %Declearing x as a variabl
P=X^5-5*X^4-6*X^3+6*X^2-3*X+1; %Equation we interest to solve
n=5; % The eq. order
Pd1 = diff(P,X,1); % first differitial of f
Pd2 = diff(P,X,2); %second differitial of f
err=0.00001; %Answear tollerance
N=100; %Max. # of Iterations
x(1)=1e-3; % Initial Value
for k=1:N
G=double(vpa(subs(Pd1,X,x(k))/subs(P,X,x(k))));
H=G^2 - double(subs(Pd2,X,x(k))) /subs(P,X,x(k));
D1= (G+sqrt((n-1)*(n*H-G^2)));
D2= (G-sqrt((n-1)*(n*H-G^2)));
D = max(D1,D2);
a=n/D;
x(k+1)=x(k)-a
Err(k) = abs(x(k+1)-x(k));
if Err(k) <=err
break
end
end输出(多项式根)：
x=
0.0010 + 0.0000i 0.1434 + 0.4661i 0.1474 + 0.4345i 0.1474 + 0.4345i 0.1474 + 0.4345i
回答 2
Stack Overflow用户
发布于 2019-12-22 17:27:51
您实际看到的是循环中产生的所有值x(k)。最后一个，0.1474 + 0.4345i是这个循环的最终结果--在给定的公差阈值中的根的近似。密码
syms X %Declaring x as a variable
P = X^5 - 5 * X^4 - 6 * X^3 + 6 * X^2 - 3 * X + 1; %Polynomial
n=5; %Degree of the polynomial
Pd1 = diff(P,X,1); %First derivative of P
Pd2 = diff(P,X,2); %Second derivative of P
err = 0.00001; %Answer tolerance
N = 100; %Maximal number of iterations
x(1) = 0; %Initial value
for k = 1:N
G = double(vpa(subs(Pd1,X,x(k)) / subs(P,X,x(k))));
H = G^2 - double(subs(Pd2,X,x(k))) / subs(P,X,x(k));
D1 = (G + sqrt((n-1) * (n * H-G^2)));
D2 = (G - sqrt((n-1) * (n * H-G^2)));
D = max(D1,D2);
a = n/D;
x(k+1) = x(k) - a;
Err(k) = abs(x(k+1)-x(k));
if Err(k) <=err
fprintf('Initial value %f, result %f%+fi', x(1), real(x(k)), imag(x(k)))
break
end
end结果：
Initial value -2.000000, result -1.649100+0.000000i如果你想得到其他的根，你必须使用其他的初始值。例如，一个人可以获得
Initial value 10.000000, result 5.862900+0.000000i
Initial value -2.000000, result -1.649100+0.000000i
Initial value 3.000000, result 0.491300+0.000000i
Initial value 0.000000, result 0.147400+0.434500i
Initial value 1.000000, result 0.147400-0.434500i这些都是多项式的零点。
当您找到另一个根时，计算下一个根的方法是通过相应的线性因子进行除法，并将循环用于生成的新多项式。请注意，这通常不是很容易处理，因为舍入错误可能会对结果产生很大影响。
Stack Overflow用户
发布于 2019-12-24 10:39:43
现有代码的问题
您没有将Laguerre方法正确地实现为复数的方法。分母候选D1,D2是一般复数，使用简单的max是不可取的，它只对实际输入有合理的结果。目的是让a=n/D成为这两个变体中较小的一个，这样我们就必须寻找绝对值更大的D in [D1,D2]。如果像C中那样有条件赋值，这将类似于
D = (abs(D_1)>abs(D2)) ? D1 : D2;由于不存在这种情况，所以必须使用具有类似结果的命令。
D = D1; if (abs(D_1)<abs(D2)) D=D2; end得到的逼近点序列是
x(0) = 0.0010000
x(1) = 0.143349512707684+0.466072958423667i
x(2) = 0.164462212064089+0.461399841949893i
x(3) = 0.164466373475316+0.461405404094130i在一个点上，我们不能期望根近似下的(残差)多项式值会大幅度减少。通过对多项式和表达式中较大项的加减，得到接近于零的数值。在这些灾难性的取消事件中丢失的准确性无法恢复。
有效为零的多项式值的阈值可以估计为double型的机器常数乘以多项式值，其中所有的系数和评价点都被它们的绝对值所取代。此测试在代码中主要用于避免零分或接近零的除法.
找到所有的根
一种方法是将该方法应用于包含所有根的圆上足够多的初始点，并在收敛太慢的情况下对早期终止有一些严格的规则。我们必须使根的列表找到唯一的，但保持多样性，.
另一种标准方法是采用通货紧缩，即把找到的根的线性因子除以。这在低度下效果很好。
不需要较慢的符号操作，因为有一些函数直接工作在系数数组上，例如polyval和polyder。利用deconv函数可以实现除以余数的通货紧缩。
对于实多项式，我们知道根的复共轭也是根。因此，用紧缩多项式初始化下一次迭代。
其他要点：
如果您不使用
- ，则创建数组是没有意义的，特别是对于
single类型而言，这是没有意义的。
示例的根
实现所有这些，我得到一个日志
x(0) = 0.001000000000000+0.000000000000000i, |Pn(x(0))| = 0.99701
x(1) = 0.143349512707684+0.466072958423667i, |dx|= 0.48733
x(2) = 0.164462212064089+0.461399841949893i, |dx|=0.021624
x(3) = 0.164466373475316+0.461405404094130i, |dx|=6.9466e-06
root found x=0.164466373475316+0.461405404094130i with value P0(x)=-2.22045e-16+9.4369e-16i
Deflation
x(0) = 0.164466373475316-0.461405404094130i, |Pn(x(0))| = 2.1211e-15
root found x=0.164466373475316-0.461405404094130i with value P0(x)=-2.22045e-16-9.4369e-16i
Deflation
x(0) = 0.164466373475316+0.461405404094130i, |Pn(x(0))| = 4.7452
x(1) = 0.586360702193454+0.016571894375927i, |dx|= 0.61308
x(2) = 0.562204173408499+0.000003168181059i, |dx|=0.029293
x(3) = 0.562204925474889+0.000000000000000i, |dx|=3.2562e-06
root found x=0.562204925474889+0.000000000000000i with value P0(x)=2.22045e-16-1.33554e-17i
Deflation
x(0) = 0.562204925474889-0.000000000000000i, |Pn(x(0))| = 7.7204
x(1) = 3.332994579372812-0.000000000000000i, |dx|= 2.7708
root found x=3.332994579372812-0.000000000000000i with value P0(x)=6.39488e-14-3.52284e-15i
Deflation
x(0) = 3.332994579372812+0.000000000000000i, |Pn(x(0))| = 5.5571
x(1) = -2.224132251798332+0.000000000000000i, |dx|= 5.5571
root found x=-2.224132251798332+0.000000000000000i with value P0(x)=-3.33067e-14+1.6178e-15i修改后的代码
P = [1, -2, -6, 6, -3, 1];
P0 = P;
deg=length(P)-1; % The eq. degree
err=1e-05; %Answer tolerance
N=10; %Max. # of Iterations
x=1e-3; % Initial Value
for n=deg:-1:1
dP = polyder(P); % first derivative of P
d2P = polyder(dP); %second derivative of P
fprintf("x(0) = %.15f%+.15fi, |Pn(x(0))| = %8.5g\n", real(x),imag(x), abs(polyval(P,x)));
for k=1:N
Px = polyval(P,x);
dPx = polyval(dP,x);
d2Px = polyval(d2P,x);
if abs(Px) < 1e-14*polyval(abs(P),abs(x))
break % if value is zero in relative accuracy
end
G = dPx/Px;
H=G^2 - d2Px / Px;
D1= (G+sqrt((n-1)*(n*H-G^2)));
D2= (G-sqrt((n-1)*(n*H-G^2)));
D = D1;
if abs(D2)>abs(D1) D=D2; end % select the larger denominator
a=n/D;
x=x-a;
fprintf("x(%d) = %.15f%+.15fi, |dx|=%8.5g\n",k,real(x),imag(x), abs(a));
if abs(a) < err*(err+abs(x))
break
end
end
y = polyval(P0,x); % check polynomial value of the original polynomial
fprintf("root found x=%.15f%+.15fi with value P0(x)=%.6g%+.6gi\n", real(x),imag(x),real(y),imag(y));
disp("Deflation");
[ P,R ] = deconv(P,[1,-x]); % division with remainder
x = conj(x); % shortcut for conjugate pairs and clustered roots
endhttps://stackoverflow.com/questions/59418490
复制相似问题
腾讯云开发者
