第10章多项式 10.1根 找出多项式的根,即多项式为零的值,可能是许多学科共同的问题,。 MATLAB求解 这个问题,并提供其它的多项式操作工具。在 MATLAB里,多项式由一个行向量表示,它 的系数是按降序排列。例如,输入多项式x4-12x3+0x2+25x+116 》p=[1-12025116 0 5116 注意,必须包括具有零系数的项。除非特别地辨认, MATLAB无法知道哪一项为零。 给出这种形式,用函数 roots找出一个多项式的根。 p) 11.7473 2.7028 1.2251+1.4672i 因为在 MATLAB中,无论是一个多项式,还是它的根,都是向量, MATLAB按惯例 规定,多项式是行向量,根是列向量。给出一个多项式的根,也可以构造相应的多项式。在 MATLAB中,命令poly执行这个任务。 》pp=poly(r Columns 1 through 4 0.1200 0.0000 0.2500 Column 5 1.1600+0.0000i >)pp=real(pp)%throw away spurious imaginary part 1.0000-12.00000.000025.0000116.0000
第 10 章 多 项 式 10.1 根 找出多项式的根，即多项式为零的值，可能是许多学科共同的问题，。MATLAB 求解 这个问题，并提供其它的多项式操作工具。在 MATLAB 里，多项式由一个行向量表示，它 的系数是按降序排列。例如，输入多项式 x 4－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
因为MA∏LAB无隙地处理复数,当用根重组多项式时,如果一些根有虚部,由于截断 误差,则poly的结果有一些小的虚部,这是很普通的。消除虚假的虚部,如上所示,只要 使用函数rea抽取实部。 10.2乘法 函数com支持多项式乘法(执行两个数组的卷积)。考虑两个多项式a(x)=x3+2x2+3x+ 4和b(x)=x3+4x2+9x+16的乘积 [1234];b=[14916] 62050758464 结果是dx)=x6+6x5+20x4+50x3+75x2+84x+64。两个以上的多项式的乘法需要重复 使用 10.3加法 对多项式加法, MATLAB不提供一个直接的函数。如果两个多项式向量大小相同,标 准的数组加法有效。把多项式a(x)与上面给出的b(x)相加 >) d=a+b 结果是dx)=2x3+6x2+12x+20。当两个多项式阶次不同,低阶的多项式必须用首零填 补,使其与高阶多项式有同样的阶次。考虑上面多项式c和d相加 e=c+[000d 84 结果是x)=x+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
因为 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 narging=mmpadd(c, -d) 16 结果是g(x)=x6+6x3+20x4+48x3+69x2+72x+44 10.4除法 在一些特殊情况,一个多项式需要除以另一个多项式。在 MATLAB中,这由函数 decoy 完成。用上面的多项式b和c 》[q,r= decon(c,b) 这个结果是b被c除,给出商多项式q和余数r,在现在情况下r是零,因为b和q的 乘积恰好是c 10.5导数
% 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为多项式求导提供了函数 polder g >>h=polder(g) 6308014413872 10.6估值 根据多项式系数的行向量,可对多项式进行加,减,乘,除和求导,也应该能对它们进 行估值。在 MATLAB中,这由函数 polyol来完成 =linspace(-1, 3) choose 100 data points between-land 3 >v=polyval(p, x) 计算x值上的p(x),把结果存在v里。然后用函数plot绘出结果。 >) plot(x, v), title(x3+4x 2-7x-10), xlabel(x) x^3+4x^2-7x-10 10 图10.1多项式估值
由于一个多项式的导数表示简单，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 ') -1 0 1 2 3 -20 -10 0 10 20 30 40 x^3+4x^2-7x-10 x 图 10.1 多项式估值
10.7有理多项式 在许多应用中,例如富里哀( Fourier),拉普拉斯( Laplace)和Z变换,出现有理多项式或 两个多项式之比。在 MATLAB中,有理多项式由它们的分子多项式和分母多项式表示。对 有理多项式进行运算的两个函数是 residue和 polder函数 residue执行部分分式展开 》num=10*12]; numerator polynomial >)den=poly(l-1;-3; -4); denominator polynomial >[res, poles, k]=residue(num, den) -6.6667 5.0000 1.6667 -4.0000 3.0000 1.0000 结果是余数、极点和部分分式展开的常数项。上面的结果说明了该问题: 10(s+2) -66667 16667 (S+1)(S+3)(s+4)S+4s+3s+ 这个函数也执行逆运算。 >[n, d]=residue(res, poles, k) 0.000010.000020.0000 1.00008.000019.000012.0000 -4.0000 3.0000 -1.0000 在截断误差内,这与我们开始时的分子和分母多项式一致。 residue也能处理重极点的 情况,尽管这里没有考虑
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 = [ ] 结果是余数、极点和部分分式展开的常数项。上面的结果说明了该问题： 10 2 1 3 4 6 6667 4 5 3 16667 1 0 ( ) ( )( )( ) s . . s s s s s s + + + + = − + + + + + + 这个函数也执行逆运算。 » [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 也能处理重极点的 情况，尽管这里没有考虑
正如前面所述,函数 polder,对多项式求导。除此之外,如果给出两个输入,则它对 有理多项式求导 >[b, a]=polder(num, den) 20-140-320-260 l16102328553456144 该结果证实: 10(s+2) 20s3-140s2-32s-260 ds(s+1)(S+3)s+4)s+16s3+102s+328s53+5532+456s+144 10.8M文件举例 本章说明了在精通 MATLAB工具箱中两个函数。这些函数说明了本章所论述的多项 式概念和如何写M文件函数。关于M文件的更多信息,参阅第8章 在讨论M文件函数的内部结构之前,我们考虑这些函数做些什么。 )n o earlier data 0.000010.000020.0000 b 20-140-320-260 amesim(n) gible leading term ans = 10.000020.0000 >)mmp2str(b) convert polynomial to string 20s^3-140s^2-320s^1-260 》mmp2str(b,x) 20x^3-140x^2-320x^1-260
正如前面所述，函数 polyder，对多项式求导。除此之外，如果给出两个输入，则它对 有理多项式求导。 » [b , a]=polyder(num , den) b = -20 -140 -320 -260 a = 1 16 102 328 553 456 144 该结果证实： d ds s s s s s s s s s s s s s { ( ) ( )( )( ) } 10 2 1 3 4 20 140 32 260 16 102 328 553 456 144 3 2 6 5 4 3 2 + + + + = − − − − + + + + + + 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) 20°(S^3+7s^2+16s^1+13) r(b,x',1) 20°(x~3+7x^2+16x^1+13) 这里函数 mmpsim删除在多项式n中近似为零的第一个系数。函数mmp2st把数值多 项式变换成等价形式的字符串表达式。该两个函数的主体是: function y o 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 o 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-XI make sure input is a row ifind(abs(x)<.99&abs(x)<tol); find insignificant indices x(i=zeros(l, length(i)) set them to zero find significant indices if isempty(i) extreme case: nothing left else 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 2+ 3s+ 4 o MMP2STR(P, V) generates the string using the variable V instead ofs. MMP2STR(2 4, z)becomes '2z2+ 3z+4 o MMP2STR(P, V, 1)factors the polynomial into the product of a
» 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
o constant and a monic polynomial %MMP2STR([234],[],1) becomes2(s^2+1.5s+2) Copyright(c)1996 by Prentice-Hall, Inc if nargin0, pm='+ else, pm=[]; end
% 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 nargin0 , pm= ' + ' ; else, pm=[ ]; end
s=Is pm pp pe]; add final terms 10.9小结 下列表格概括了在本章所讨论的多项式操作特性。 表10.1 项式函数 conv(a, b) 乘法 Ia, rI=decon(a, b) 除法 poly(r 用根构造多项式 polder(a) 对多项式或有理多项式求导 polyfit(x, y, n) 多项式数据拟合 polyval(p, x) 计算x点中多项式值 Ir, p, k]=residue(a, b) 部分分式展开式 [a, b] 部分分式组合 求多项式的根 表10.2 精通 MATLAB多项式操作 多项式向量到字符串变换,a(s) 多项式向量到字符串变换,a(x) mmp2str(a,x, 1) 常数和符号多项式变换 多项式加法 mmpsimt a 多项式简化
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) 多项式简化