木泥的BLOG

Chapter 4: Numerical computation of MATLAB

 

一、多项式(Polynomial)`

1．多项式的表达与创建(Expression and Creating of polynomial)

(1) 多项式的表达（expression of polynomial）_

Matlab用行矢量表达多项式系数(Coefficient)，各元素按变量的降幂顺序排列，如多项式为：

P(x)=a₀xⁿ+a₁x^n-1+a₂x^n-2…a_n-1x+a_n

则其系数矢量(Vector of coefficient)为：P=[a₀ a₁ … a_n-1 a_n]

如将根矢量(Vector of root)表示为：

ar=[ ar₁ ar₂ … ar_n]

则根矢量与系数矢量之间关系为：

（x-ar₁）(x- ar₂) … (x- ar_n)= a₀xⁿ+a₁x^n-1+a₂x^n-2…a_n-1x+a_n

(2)多项式的创建（polynomial creating）

a)系数矢量的直接输入法

利用poly2sym函数直接输入多项式的系数矢量，就可方便的建立符号形式的多项式。

例：创建多项式x³-4x²+3x+2

poly2sym([1 -4 3 2])

ans =

x^3-4*x^2+3*x+2

POLY Convert roots to polynomial.

    POLY(A), when A is an N by N matrix, is a row vector with

    N+1 elements which are the coefficients of the characteristic polynomial, DET(lambda*EYE(SIZE(A)) - A) .

 

POLY(V), when V is a vector, is a vector whose elements are

    the coefficients of the polynomial whose roots are the

    elements of V . For vectors, ROOTS and POLY are inverse

    functions of each other, up to ordering, scaling, and

    roundoff error.

b) 由根矢量创建多项式

通过调用函数 p=poly(ar)产生多项式的系数矢量, 再利用poly2sym函数就可方便的建立符号形式的多项式。

 

注：（1）根矢量元素为n ，则多项式系数矢量元素为n+1；

（2）函数poly2sym(pa) 把多项式系数矢量表达成符号形式的多项式，缺省情况下自变量符号为x，可以指定自变量。

（3）使用简单绘图函数ezplot可以直接绘制符号形式多项式的曲线。

 

例 1：由根矢量创建多项式。将多项式(x-6)(x-3)(x-8)表示为系数形式

 a=[6  3  8]                    %根矢量

pa=poly(a)                       %求系数矢量

ppa=poly2sym(pa)                %以符号形式表示原多项式

ezplot(ppa,[-50,50])

pa =

    1   -17    90  -144

ppa =

x^3-17*x^2+90*x-144

注：含复数根的根矢量所创建的多项式要注意：

（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；

   （2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。

进行多项式的求根运算时，有两种方法，一是直接调用求根函数roots，poly和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。

例 3： 由给定复数根矢量求多项式系数矢量。

r=[-0.5  -0.3+0.4i  -0.3-0.4i];

p=poly(r)

pr=real(p)

ppr=poly2sym(pr)

p =

    1.0000    1.1000    0.5500    0.1250

pr =

    1.0000    1.1000    0.5500    0.1250

ppr =

x^3+11/10*x^2+11/20*x+1/8

 

c) 特征多项式输入法

用poly函数可实现由矩阵的特征多项式系数创建多项式。

条件：特征多项式系数矢量的第一个元素必须为一。

例 2： 求三阶方阵A的特征多项式系数，并转换为多项式形式。

a=[6  3  8;7  5  6; 1  3  5]

Pa=poly(a)       %求矩阵的特征多项式系数矢量

Ppa=poly2sym(pa)

Pa =

   1.0000  -16.0000   38.0000  -83.0000

Ppa =

x^3-17*x^2+90*x-144

 

注：n 阶方阵的特征多项式系数矢量一定是n +1阶的。

 

注：（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；

   （2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。

进行多项式的求根运算时，有两种方法，一是直接调用求根函数roots，poly和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。

例 4： 将多项式的系数表示形式转换为根表现形式。

求 x³-6x²-72x-27的根

a=[1  -6  -72  -27]

r=roots(a)

r =

   12.1229

   -5.7345

   -0.3884

MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。

>> 

1. 多项式的乘除运算(Multiplication and division of polynomial)

多项式乘法用函数conv(a,b)实现， 除法用函数deconv(a,b)实现。

例1：a(s)=s²+2s+3,  b(s)=4s²+5s+6,计算 a(s)与 b(s)的乘积。

a=[1 2 3];  b=[4 5 6];

c=conv(a,b)

cs=poly2sym(c,’s’)

c =

     4    13    28    27    18

cs =

4*s^4+13*s^3+28*s^2+27*s+18

例2： 展开(s²+2s+2)(s+4)(s+1)  （多个多项式相乘）

c=conv([1,2,2],conv([1,4],[1,1]))

cs=poly2sym(c,’s’)            %（指定变量为s）

c =

     1     7    16    18     8

cs =

s^4+7*s^3+16*s^2+18*s+8

 

例2：求多项式s^4+7*s^3+16*s^2+18*s+8分别被(s+4),(s+3)除后的结果。

c=[1  7  16  18  8];

[q1,r1]=deconv(c,[1,4])      %q—商矢量，  r—余数矢量

[q2,r2]=deconv(c,[1,3])

cc=conv(q2,[1,3])            %对除(s+3)结果检验

test=((c-r2)==cc)

q1 =

   1     3     4     2

r1 =

    0     0     0     0     0

q2 =

    1     4     4     6

r2 =

    0     0     0     0   -10

cc =

    1     7    16    18    18

test =

    1     1     1     1     1

1. 其他常用的多项式运算命令(Other computation command of polynomial)

pa=polyval(p,s)    按数组运算规则计算给定s时多项式p的值。

pm=polyvalm(p,s)  按矩阵运算规则计算给定s时多项式p的值。

[r,p,k]=residue(b,a)  部分分式展开，b,a分别是分子分母多项式系数矢量，r,p,k分别是留数、极点和直项矢量

p=polyfit(x,y,n)     用n阶多项式拟合x，y矢量给定的数据。

polyder(p)         多项式微分。

注： 对于多项式b(s)与不重根的n阶多项式a(s)之比，其部分分式展开为：  

式中：p₁,p₂,…,p_n称为极点(poles)，r₁,r₂,…,r_n 称为留数(residues)，k(s)称为直项(direct terms)，假如a(s)含有m重根p_j，则相应部分应写成：

RESIDUE Partial-fraction expansion (residues).

    [R,P,K] = RESIDUE(B,A) finds the residues, poles and direct term of a partial fraction expansion of the ratio of two polynomials B(s)/A(s). If there are no multiple roots,

B(s)    R(1)      R(2)         R(n)

        ----  = --------  +  -------- + ... +  --------  +  K(s)

        A(s)   s - P(1)    s - P(2)       s - P(n)

Vectors B and A specify the coefficients of the numerator and denominator polynomials in descending powers of s.  The residues

are returned in the column vector R, the pole locations in column vector P, and the direct terms in row vector K. The number of poles is n = length(A)-1 = length(R) = length(P). The direct term  coefficient vector is empty if length(B) < length(A), otherwise

    length(K) = length(B)-length(A)+1.

If P(j) = ... = P(j+m-1) is a pole of multplicity m, then the expansion includes terms of the form

 R(j)     R(j+1)         R(j+m-1)

-------- + ------------   + ... + ------------

s - P(j)   (s - P(j))^2     (s - P(j))^m

[B,A] = RESIDUE(R,P,K), with 3 input arguments and 2 output arguments, converts the partial fraction expansion back to the polynomials with coefficients in B and A.

 

例3：对 (3x⁴+2x³+5x²+4x+6)/(x⁵+3x⁴+4x³+2x²+7x+2)  做部分分式展开

a=[1  3  4  2  7  2];

b=[3  2  5  4  6];

[r,s,k]=residue(b,a)

r =

   1.1274 + 1.1513i

   1.1274 - 1.1513i

  -0.0232 - 0.0722i

  -0.0232 + 0.0722i

   0.7916         

s =

  -1.7680 + 1.2673i

  -1.7680 - 1.2673i

   0.4176 + 1.1130i

   0.4176 - 1.1130i

  -0.2991         

k =

[]   (分母阶数高于分子阶数时，k将是空矩阵，表示无此项)

 

例5：对一组实验数据进行多项式最小二乘拟合(least square fit)

x=[1  2  3  4  5];                  % 实验数据

y=[5.5  43.1  128  290.7  498.4];

p=polyfit(x,y,3)                    %做三阶多项式拟合

x2=1:.1:5;

y2=polyval(p,x2);         % 根据给定值计算多项式结果

plot(x,y,’o’,x2,y2)