【 * 例 3.8.2 -1 】求 的“商”及“余”多项式。
p1=conv([1,0,2],conv([1,4],[1,1])); % 计算分子多项式
p2=[1 0 1 1]; % 注意缺项补零
[q,r]=deconv(p1,p2);
cq=' 商多项式为 '; cr=' 余多项式为 ';
disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')])
商多项式为 s + 5
余多项式为 5 s^2 + 4 s + 3
【 * 例 3.8.2 -2 】两种多项式求值指令的差别。
S=pascal(4) % 生成一个 4 阶方阵
P=poly(S);PP=poly2str(P,'s')
PA=polyval(P,S) % 独立变量取数组 S 元素时的多项式值
PM=polyvalm(P,S) % 独立变量取矩阵 S 时的多项式值
S =
1 1 1 1
1 2 3 4
1 3 6 10
1 4 10 20
PP =
s^4 - 29 s^3 + 72 s^2 - 29 s + 1
PA =
1.0e+004 *
0.0016 0.0016 0.0016 0.0016
0.0016 0.0015 -0.0140 -0.0563
0.0016 -0.0140 -0.2549 -1.2089
0.0016 -0.0563 -1.2089 -4.3779
PM =
1.0e-011 *
-0.0077 0.0053 -0.0096 0.0430
-0.0068 0.0481 -0.0110 0.1222
0.0075 0.1400 -0.0095 0.2608
0.0430 0.2920 -0.0007 0.4737
【 * 例 3.8.2 -3 】部分分式展开。
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 =
[]
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 =
162050758464
结果是c(x)=x6+6x5+20x4+50x3+75x2+84x+64。两个以上的多项式的乘法需要重复使用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)<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=;%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) | 多项式简化 |
文章评论(0条评论)
登录后参与讨论