matlab多项式

10.1 根

找出多项式的根，即多项式为零的值，可能是许多学科共同的问题。MATLAB求解这个问题，并提供其它的多项式操作工具。在

MATLAB里，多项式由一个行向量表示，它的系数是按降序排列。例如，输入多项式x4－12x3＋0x2＋25x＋116

» p=[1-12025116]

p =

1-12025116

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

。

» r=roots(p)

r =

11.7473

2.7028

-1.2251 + 1.4672i

-1.2251 - 1.4672i

因为在MATLAB中，无论是一个多项式，还是它的根，都是向量，MATLAB按惯例规定，多项式是行向量，根是列向量。给出一

» pp=poly(r)

pp =

1.0e+002 *

Columns 1 through 4 

0.0100-0.12000.00000.2500

Column 5 

1.1600 + 0.0000i 

» pp=real(pp) %throw away spurious imaginary part

pp =

1.0000-12.00000.000025.0000116.0000

因为MATLAB无隙地处理复数，当用根重组多项式时，如果一些根有虚部，由于截断误差，则poly的结果有一些小的虚部，这是很普

通的。消除虚假的虚部，如上所示，只要使用函数real抽取实部。

10.2 乘法

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

» a=[1234] ;b=[14916];

» c=conv(a , b)

c =

1620507584

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

10.3 加法

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

b(x)相加。

» d=a+b

d =

261220

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

上面多项式c和d相加：

» e=c+[000d]

e =

162052819684

结果是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 =

162052819684

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

»g=mmpadd(c , -d)

g =

162048697244

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

10.4 除法

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

» [q , r]=deconv(c , b)

q =

1234

r =

0000000

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

10.5 导数

由于一个多项式的导数表示简单，MATLAB为多项式求导提供了函数polyder。

» g

g =

162048697244

» h=polyder(g)

h =

6308014413872

10.6估值

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

成。

» x=linspace(-1, 3) ; %choose 100 data points between -1and 3.

» p=[14-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*[12] ; %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.000010.000020.0000

d =

1.00008.000019.000012.0000

» roots(d)

ans =

-4.0000

-3.0000

-1.0000

在截断误差内，这与我们开始时的分子和分母多项式一致。residue也能处理重极点的情况，尽管这里没有考虑。

正如前面所述，函数polyder，对多项式求导。除此之外，如果给出两个输入，则它对有理多项式求导。

» [b , a]=polyder(num , den)

b =

-20-140-320-260

a =

116102328553456144

该结果证实：

10.8 M文件举例

本章说明了在精通MATLAB工具箱 中两个函数。这些函数说明了本章所论述的多项式概念和如何写M文件函数。关于M文件的更多信

息，参阅第8章。

在讨论M文件函数的内部结构之前，我们考虑这些函数做些什么。

» n%earlier data

n =

0.000010.000020.0000

» b%earlier data

b =

-20-140-320-260

» mmpsim(n)%strip away negligible leading term

ans =

10.000020.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)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 = [234] becomes the string ' 2s^2 + 3s + 4 '

%

%MMP2STR(P,V) generates the string using the variable V

%instead of s. MMP2STR([234],' z ') becomes ' 2z^2 + 3z + 4 '

%

%MMP2STR(P,V,1) factors the polynomial into the product of a

%constant and a monic polynomial.

%MMP2STR([234],[],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=[spmppv' ^ ' sprintf(' %.0f ',n-i)] ;end

end

if p(n)~ =0,pp=sprintf(' %.4g ',abs(p(n))); else,pp=[];end

if p(n)<0 , pm= ' -' ; elseifp(n)>0 , pm= ' +' ; else,pm=[];end

s=[spmpppe];%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)

多项式简化