|
10.1 根 找出多項式的根，即多項式為零的值，可能是許多學(xué)科共同的問題，。MATLAB求解這個問題，并提供其它的多項式操作工具。在MATLAB里，多項式由一個行向量表示，它的系數(shù)是按降序排列。例如，輸入多項式x4－12x3＋0x2＋25x＋116 p=[1-12025116] p = 1-12025116 注意，必須包括具有零系數(shù)的項。除非特別地辨認(rèn)，MATLAB無法知道哪一項為零。給出這種形式，用函數(shù)roots找出一個多項式的根。 r=roots(p) r = 11.7473 2.7028 -1.2251 + 1.4672i -1.2251 - 1.4672i 因為在MATLAB中，無論是一個多項式，還是它的根，都是向量，MATLAB按慣例規(guī)定，多項式是行向量，根是列向量。給出一 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無隙地處理復(fù)數(shù)，當(dāng)用根重組多項式時，如果一些根有虛部，由于截斷誤差，則poly的結(jié)果有一些小的虛部，這是很普通的。消除虛假的虛部，如上所示，只要使用函數(shù)real抽取實部。 10.2 乘法 函數(shù)conv支持多項式乘法(執(zhí)行兩個數(shù)組的卷積)。考慮兩個多項式a(x)=x3＋2x2＋3x＋4和b(x)= x3＋4x2＋9x＋16的乘積： a=[1234] ;b=[14916]; c=conv(a , b) c = 162050758464 結(jié)果是c(x)=x6＋6x5＋20x4＋50x3＋75x2＋84x＋64。兩個以上的多項式的乘法需要重復(fù)使用conv。 10.3 加法 對多項式加法，MATLAB不提供一個直接的函數(shù)。如果兩個多項式向量大小相同，標(biāo)準(zhǔn)的數(shù)組加法有效。把多項式a(x)與上面給出的b(x)相加。 d=a+b d = 261220 結(jié)果是d(x)= 2x3＋6x2＋12x＋20。當(dāng)兩個多項式階次不同，低階的多項式必須用首零填補(bǔ)，使其與高階多項式有同樣的階次?？紤]上面多項式c和d相加： e=c+[000d] e = 162052819684 結(jié)果是e(x)= x6＋6x5＋20x4＋52x3＋81x2＋96x＋84。要求首零而不是尾零，是因為相關(guān)的系數(shù)象x冪次一樣，必須整齊。 如果需要，可用一個文件編輯器創(chuàng)建一個函數(shù)M文件來執(zhí)行一般的多項式加法。精通MATLAB工具箱包含下列實現(xiàn)： 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 現(xiàn)在，為了闡述mmpadd的使用，再考慮前一頁的例子。 f=mmpadd(c,d) f = 162052819684 它與上面的e相同。當(dāng)然，mmpadd也用于減法。 g=mmpadd(c , -d) g = 162048697244 10.4 除法 在一些特殊情況，一個多項式需要除以另一個多項式。在MATLAB中，這由函數(shù)deconv完成。用上面的多項式b和c [q , r]=deconv(c , b) q = 1234 r = 0000000 這個結(jié)果是b被c除，給出商多項式q和余數(shù)r，在現(xiàn)在情況下r是零，因為b和q的乘積恰好是c。 10.5 導(dǎo)數(shù) 由于一個多項式的導(dǎo)數(shù)表示簡單，MATLAB為多項式求導(dǎo)提供了函數(shù)polyder。 g g = 162048697244 h=polyder(g) h = 6308014413872 10.6估值 根據(jù)多項式系數(shù)的行向量，可對多項式進(jìn)行加，減，乘，除和求導(dǎo)，也應(yīng)該能對它們進(jìn)行估值。在MATLAB中，這由函數(shù)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)，把結(jié)果存在v里。然后用函數(shù)plot繪出結(jié)果。 plot(x , v),title(' x^3+4x^2-7x-10 '),xlabel(' x ')
圖10.1多項式估值 10.7有理多項式 在許多應(yīng)用中，例如富里哀(Fourier)，拉普拉斯(Laplace)和Z變換，出現(xiàn)有理多項式或兩個多項式之比。在MATLAB中，有理多項式由它們的分子多項式和分母多項式表示。對有理多項式進(jìn)行運(yùn)算的兩個函數(shù)是residue和polyder。函數(shù)residue執(zhí)行部分分式展開。 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 = [] 結(jié)果是余數(shù)、極點(diǎn)和部分分式展開的常數(shù)項。上面的結(jié)果說明了該問題：
這個函數(shù)也執(zhí)行逆運(yùn)算。 [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 在截斷誤差內(nèi)，這與我們開始時的分子和分母多項式一致。residue也能處理重極點(diǎn)的情況，盡管這里沒有考慮。 正如前面所述，函數(shù)polyder，對多項式求導(dǎo)。除此之外，如果給出兩個輸入，則它對有理多項式求導(dǎo)。 [b , a]=polyder(num , den) b = -20-140-320-260 a = 116102328553456144 該結(jié)果證實：
10.8 M文件舉例 本章說明了在精通MATLAB工具箱 中兩個函數(shù)。這些函數(shù)說明了本章所論述的多項式概念和如何寫M文件函數(shù)。關(guān)于M文件的更多信息，參閱第8章。 在討論M文件函數(shù)的內(nèi)部結(jié)構(gòu)之前，我們考慮這些函數(shù)做些什么。 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 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) 這里函數(shù)mmpsim刪除在多項式n中近似為零的第一個系數(shù)。函數(shù)mmp2str把數(shù)值多項式變換成等價形式的字符串表達(dá)式。該兩個函數(shù)的主體是： 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 = [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 小結(jié) 下列表格概括了在本章所討論的多項式操作特性。 表10.1
表10.2
|
||||||||||||||||||||||||||||||||
|
|
