c4_6.m

这是我从一本基于MATLAB实现的数值分析算法的书中自带的光盘程序

%C4_6.m
%用劈因子法求多项式方程的一对共轭复根
function C4_6
n=input('The order of polynomial?');
a(1)=input('a(0)=');
a(2)=input('a(1)=');
a(3)=input('a(2)=');
a(4)=input('a(3)=');
a(5)=input('a(4)=');
ep=input('Tolerance:ep?');
p0=input('p0=');
q0=input('q0=');
b(1)=99;
b(2)=99;
while abs(b(1))>ep |abs(b(2))>ep
%求简化的多项式的系数b
      b(n+1)=a(n+1);
      b(n)=a(n)-p0*b(n+1);
      for i=n-1:-1:2
          b(i)=a(i)-p0*b(i+1)-q0*b(i+2);
      end
      b(1)=a(1)-q0*b(3);
%对p的导数bp
      bp(n+1)=0;
      bp(n)=-b(n+1)-p0*bp(n+1);
      for i=n-1:-1:2
          bp(i)=-b(i+1)-p0*bp(i+1)-q0*bp(i+2);
      end
      bp(1)=-q0*bp(3);
%对q的导数bq
      bq(n+1)=0;
      bq(n)=0;
      for i=n-1:-1:2
          bq(i)=-b(i+2)-p0*bq(i+1)-q0*bq(i+2);
      end
      bq(1)=-b(3)-q0*bq(3);
%求解线性方程组
      A=[bp(1),bq(1);bp(2),bq(2)];
      d=[-b(1),-b(2)]';
      x=A\d;
      dp=x(1);
      dq=x(2);
      p0=p0+dp;
      q0=q0+dq;
end
p=p0;q=q0;
fprintf('p=%3.1f,q=%3.1f \n',p,q);
fprintf('Quadratic coefficient=x^2+(%3.1f x)+(%3.1f) \n ',p,q );
t1(1)=-p/2;t1(2)=sqrt(-p^2+4*q)/2;
t2(1)=-p/2;t2(2)=-sqrt(-p^2+4*q)/2;
fprintf('Roots of the quadratic factor: \n');
fprintf('%3.2f+%8.6fi \n',t1(1),t1(2));
fprintf('%3.2f+%8.6fi \n',t2(1),t2(2));
fprintf('Coefficients of deflated polynomial: \n');
fprintf('order , coefficients  \n');
fprintf('%3.0f , %6.5f \n',n-4,b(n-1));
fprintf('%3.0f , %6.5f \n',n-3,b(n));
fprintf('%3.0f , %6.5f \n',n-2,b(n+1));