Matlab 版 (精华区)

发信人: bage (最近比较烦), 信区: Matlab
标  题: 精通Matlab(六)
发信站: 哈工大紫丁香 (Sun Feb  4 14:03:49 2001), 转信

发信人: Security (淼水), 信区: MathTools       
发信站: BBS 水木清华站 (Tue Jun  1 22:36:08 1999)


第10章  多 项 式


10.1  根

    找出多项式的根，即多项式为零的值，可能是许多学科共同的问题，。MATLAB求解这个问题，并提供其它的多项式操作工具。在MATLAB里，多项式由一个行向量表示，它的系数是按降序排列。例如，输入多项式x4－12x3＋0x2＋25x＋116

                   ? p=[1  -12  0  25  116]
                p =
                     1   -12     0    25   116

    注意，必须包括具有零系数的项。除非特别地辨认，MATLAB无法知道哪一项为零。给出这种形式，用函数roots找出一个多项式的根。

                        ? r=roots(p)
                        r =
                          11.7473          
                           2.7028          
                          -1.2251 + 1.4672i
                          -1.2251 - 1.4672i

    因为在MATLAB中，无论是一个多项式，还是它的根，都是向量，MATLAB按惯例规定，多项式是行向量，根是列向量。给出一个多项式的根，也可以构造相应的多项式。在MATLAB中，命令poly执行这个任务。

        ? pp=poly(r)
        pp =
          1.0e+002 *
          Columns 1 through 4 
           0.0100            -0.1200             0.0000             0.2500          
          Column 5 
           1.1600 + 0.0000i 

        ? pp=real(pp) %throw away spurious imaginary part
        pp =
            1.0000  -12.0000    0.0000   25.0000  116.0000

    因为MATLAB无隙地处理复数，当用根重组多项式时，如果一些根有虚部，由于截断误差，则poly的结果有一些小的虚部，这是很普通的。消除虚假的虚部，如上所示，只要使用函数real抽取实部。

10.2  乘法

    函数conv支持多项式乘法(执行两个数组的卷积)。考虑两个多项式a(x)=x3＋2x2＋3x＋4和b(x)= x3＋4x2＋9x＋16的乘积：

                ? a=[1  2  3  4] ;  b=[1  4  9  16];
                ? c=conv(a , b)
                c =
                     1     6    20    50    75    84    64

    结果是c(x)=x6＋6x5＋20x4＋50x3＋75x2＋84x＋64。两个以上的多项式的乘法需要重复使用conv。

10.3  加法

    对多项式加法，MATLAB不提供一个直接的函数。如果两个多项式向量大小相同，标准的数组加法有效。把多项式a(x)与上面给出的b(x)相加。

                        ? d=a+b
                        d =
                             2     6    12    20

    结果是d(x)= 2x3＋6x2＋12x＋20。当两个多项式阶次不同，低阶的多项式必须用首零填补，使其与高阶多项式有同样的阶次。考虑上面多项式c和d相加：

                ? e=c+[0  0  0  d]
                e =
                     1     6    20    52    81    96    84

    结果是e(x)= x6＋6x5＋20x4＋52x3＋81x2＋96x＋84。要求首零而不是尾零，是因为相关的系数象x幂次一样，必须整齐。
    如果需要，可用一个文件编辑器创建一个函数M文件来执行一般的多项式加法。精通MATLAB工具箱包含下列实现：

        function p=mmpadd(a,b)
        %  MMPADD Polynomial addition.
        %  MMPADD(A,B) adds the polynomial A and B

        %  Copyright (c) 1996 by Prentice Hall,Inc.

        if nargin<2
             error(' Not enough input arguments ')
        end
        a=a(:).' ;          %  make sure inputs are polynomial row vectors
        b=b(:).' ;
        na=length(a) ;       %  find lengths of a and b
        nb=length(b) ;
        p=[zeros(1,nb-na) a]+[zeros(1,na-nb) b] ;     %  add zeros as necessary

    现在，为了阐述mmpadd的使用，再考虑前一页的例子。

                ? f=mmpadd(c,d)
                f =
                     1     6    20    52    81    96    84

    它与上面的e相同。当然，mmpadd也用于减法。

                ?g=mmpadd(c , -d)
                g =
                     1     6    20    48    69    72    44

    结果是g(x)= x6＋6x5＋20x4＋48x3＋69x2＋72x＋44。

10.4  除法

    在一些特殊情况，一个多项式需要除以另一个多项式。在MATLAB中，这由函数deconv完成。用上面的多项式b和c 

                ? [q , r]=deconv(c , b)
                q =
                     1     2     3     4
                r =
                     0     0     0     0     0     0     0

    这个结果是b被c除，给出商多项式q和余数r，在现在情况下r是零，因为b和q的乘积恰好是c。

10.5  导数

    由于一个多项式的导数表示简单，MATLAB为多项式求导提供了函数polyder。
                ? g
                g =
                     1     6    20    48    69    72    44
                ? h=polyder(g)
                h =
                     6    30    80   144   138    72

10.6  估值

    根据多项式系数的行向量，可对多项式进行加，减，乘，除和求导，也应该能对它们进行估值。在MATLAB中，这由函数polyval来完成。

        ? x=linspace(-1, 3) ;      %  choose 100 data points between -1and 3.
        ? p=[1  4  -7  -10] ;    %  uses polynomial p(x) = x3＋4x2－7x－10
        ? v=polyval(p , x) ;

    计算x值上的p(x)，把结果存在v里。然后用函数plot绘出结果。

        ? plot(x , v),title(' x^3+4x^2-7x-10 '),  xlabel(' x ')


 
图10.1  多项式估值

10.7  有理多项式

    在许多应用中，例如富里哀(Fourier)，拉普拉斯(Laplace)和Z变换，出现有理多项式或两个多项式之比。在MATLAB中，有理多项式由它们的分子多项式和分母多项式表示。对有理多项式进行运算的两个函数是residue和polyder。函数residue执行部分分式展开。

                ? num=10*[1  2] ;           %  numerator polynomial
                ? den=poly([-1;  -3;  -4]) ;   %  denominator polynomial
                ? [res, poles, k]=residue(num, den)
                res =
                   -6.6667
                    5.0000
                    1.6667
                poles =
                   -4.0000
                   -3.0000
                   -1.0000
                k =
                     [  ]

    结果是余数、极点和部分分式展开的常数项。上面的结果说明了该问题：

     

这个函数也执行逆运算。

                ? [n, d]=residue(res, poles, k)
                n =
                    0.0000   10.0000   20.0000
                d =
                    1.0000    8.0000   19.0000   12.0000

                ? roots(d)
                ans =
                   -4.0000
                   -3.0000
                   -1.0000

    在截断误差内，这与我们开始时的分子和分母多项式一致。residue也能处理重极点的情况，尽管这里没有考虑。
    正如前面所述，函数polyder，对多项式求导。除此之外，如果给出两个输入，则它对有理多项式求导。

                ? [b , a]=polyder(num , den)
                b =
                   -20  -140  -320  -260
                a =
                     1    16   102   328   553   456   144

    该结果证实：

     

10.8  M文件举例

    本章说明了在精通MATLAB工具箱 中两个函数。这些函数说明了本章所论述的多项式概念和如何写M文件函数。关于M文件的更多信息，参阅第8章。
    在讨论M文件函数的内部结构之前，我们考虑这些函数做些什么。

                ? n      %  earlier data
                n =
                    0.0000   10.0000   20.0000

                ? b      %  earlier data
                b =
                   -20  -140  -320  -260

                ? mmpsim(n)     %  strip away negligible leading term
                ans =
                   10.0000   20.0000

                ? mmp2str(b)     %  convert polynomial to string
                ans =
                -20s^3 - 140s^2 - 320s^1 - 260

                ? mmp2str(b , ' x ')  
                ans =
                -20x^3 - 140x^2 - 320x^1 - 260

                ? mmp2str(b , [  ] , 1)  
                ans =
                -20*(s^3 + 7s^2 + 16s^1 + 13)

                ? mmp2str(b , ' x ' , 1)  
                ans =
                -20*(x^3 + 7x^2 + 16x^1 + 13)

    这里函数mmpsim删除在多项式n中近似为零的第一个系数。函数mmp2str把数值多项式变换成等价形式的字符串表达式。该两个函数的主体是：

        function y=mmpsim(x,tol)
        %  MMPSIM Polynomial Simplification,Strip Leading Zero Terms.
        %  MMPSIM(A) Delets leading zeros and small coefficients in the 
        %  polynomial A(s). Coefficients are considered small if their 
        %  magnitude is less than both one and norm(A)*1000*eps.
        %  MMPSIM(A,TOL) uses TOL for its smallness tolerance.

        %  Copyright (c) 1996 by Prentice-Hall,Inc.

        if nargin<2, tol=norm(x)*1000*eps; end
        x=x(:).' ;                       %  make sure input is a row 
        i=find(abs(x)<.99&abs(x)<tol) ;    %  find insignificant indices
        x(i)=zeros(1, length(i)) ;          %  set them to zero
        i=find(x~=0) ;                  %  find significant indices
        if isempty(i)
             y=0 ;                     %  extreme case: nothing left!
        else
            y=x(i(1) : length(x)) ;         %  start with first term
        end                           %  and go to end of polynomial


        function s=mmp2str(p,v,ff)
        %  MMP2STR Polynomial Vector to String Conversion.
        %  MMP2STR(P) converts the polynomial vector P into a string.
        %  For example: P = [2  3  4] becomes the string ' 2s^2 + 3s + 4 '
        %
        %  MMP2STR(P,V) generates the string using the variable V
        %  instead of s. MMP2STR([2  3  4],' z ') becomes ' 2z^2 + 3z + 4 '
        %
        %  MMP2STR(P,V,1) factors the polynomial into the product of a
        %  constant and a monic polynomial.
        %  MMP2STR([2  3  4],[  ],1) becomes ' 2(s^2 + 1.5s + 2) '

        %  Copyright (c) 1996 by Prentice-Hall,Inc.

        if nargin<3, ff=0; end            %  factored form is False
        if nargin <2, v=' s ' ; end          %  default variable is ' s '
        if isempty(v), v=' s ' ; end         %  default variable is ' s '
        v=v(1) ;                       %  variable must be scalar 
        p=mmpsim(p) ;                 %  strip insignificant terms
        n=length(p) ;
        if ff                           %  put in factored form
            K=p(1) ; Ka=abs(K) ; p=p/K;
            if abs(K-1)<1e-4
               pp=[  ]; pe=[  ] ;
            elseif abs(K+1)<1e-4
               pp=' -(' ; pe= ') ' ;
            elseif abs(Ka-round(Ka))<=1e-5*Ka
               pp=[sprintf(' %.0f ', K) '*( ' ] ; pe= ') ' ;
            else
               pp=[sprintf(' %.4g ' , K) '*(' ] ; pe= ') ' ;
            end
        else                           %  not factored form
            K=p(1);
            pp=sprintf(' %.4g ' , K) ;
            pe=[  ];
        end
        if n==1                        %  polynomial is just a constant
            s=sprintf(' %.4g ',K);
            return
        end
        s=[pp v ' ^ ' sprintf(' %.0f ',n-1)];    %  begin string construction

        for i=2:n-1                      %  catenate center terms in polynomial
           if p(i)<0, pm= ' -  ' ;  else,  if p(i)<0,pm= ' ;  end
           if p(i)= =1,pp=[  ] ; else,  pp=sprintf(' %.4g ', abs(p(i))) ;  end
           if p(i)~ =0,s=[s  pm  pp  v  ' ^ ' sprintf(' %.0f ',n-i)] ;  end
        end

        if p(n)~ =0,pp=sprintf(' %.4g ',abs(p(n))); else,  pp=[  ];  end
        if p(n)<0 , pm= ' -  ' ; elseif  p(n)>0 , pm= ' +  ' ; else,  pm=[  ];  end
        s=[s  pm  pp  pe];             %  add final terms


10.9  小结

    下列表格概括了在本章所讨论的多项式操作特性。

表10.1
多   项   式   函   数
conv(a, b)      乘法
[q, r]=deconv(a, b)     除法
poly(r) 用根构造多项式
polyder(a)      对多项式或有理多项式求导
polyfit(x, y, n)        多项式数据拟合
polyval(p, x)   计算x点中多项式值
[r, p, k]=residue(a, b) 部分分式展开式
[a, b]=residue(r, p, k) 部分分式组合
roots(a)        求多项式的根


表10.2
精 通 MATLAB 多 项 式 操 作
mmp2str(a)      多项式向量到字符串变换，a(s)
mmp2str(a, ' x ')       多项式向量到字符串变换，a(x)
mmp2str(a, ' x ', 1)    常数和符号多项式变换
mmpadd(a, b)    多项式加法
mmpsim(a)       多项式简化



--

--
☆ 来源:．哈工大紫丁香 bbs.hit.edu.cn．[FROM: bage.bbs@smth.org]