Matlab 版 (精华区)

发信人: xuedy (dingyu), 信区: Matlab
标  题: Re: 我想知道roots(),poly()函数的算法我该怎么...
发信站: 哈工大紫丁香 (2002年03月11日18:37:14 星期一), 站内信件


【 在 summersun (羊咩咩) 的大作中提到: 】
: 【 在 Christy (绿叶~~捣鼓六仙捣毁仙) 的大作中提到: 】
: : help sort
: : %就可以把这个函数的函数体打印出来了
: : 【 在 summersun (羊咩咩) 的大作中提到: 】
: 
: thank you
: 

其实可以看出，roots和poly实际上是基于内在函数eig的，所以你如果想在底层使用，
则应该将LaPack中的相应源程序也考虑在内。
poly函数写的不好，精度不高，我在科学运算书中给出了一个基于Fadeev算法的MATLAB
函数，效果好，速度也快。
你可以自己比较一下：
例子，
>> v=[1 2 3 4 5]; V=vander(v);
   v1=poly(V); v2=poly1(V) % poly1是我编的函数
   v3=poly(sym(V)) %符号运算工具箱解
现在验证Hamilton-Cailey定理：
>> va1=polyvalm(v1,V)
   va2=polyvalm(v2,V)
清单如下，其好处还在于它不调用任何内在函数：
function c=poly1(A)
%POLY1 函数用 Fadeev-Faddeva 算法来取矩阵的特征多项式系数。
%   它完全可以取代 MATLAB 本身提供的 POLY 函数，且运算精度更高。
%Designed by Prof D Xue (c) 2000
[nr,nc]=size(A);
if nc==nr % 给出若为方阵，则用Fadeev算法求特征多项式
   I=eye(nc); R=I; c=[1 zeros(1,nc)];
   for k=1:nc
      c(k+1)=-1/k*trace(A*R);
      R=A*R+c(k+1)*I;
   end
elseif (nr==1 | nc==1) % 给出为向量时，构造矩阵
   A=A(isfinite(A)); n=length(A); % 除去非数或无界的特征根
   c=[1 zeros(1,n)];
   for j=1:n
      c(2:(j+1))=c(2:(j+1))-A(j).*c(1:j);
   end
else   % 参数有误则给出错误信息
   error('Argument must be a vector or a square matrix.')
end


--
※ 来源:·哈工大紫丁香 bbs.hit.edu.cn·[FROM: 218.25.30.141]