第四章 MATLAB的数值计算功能 Chapter 4: Numerical computation of MATLaB 、多项式( Polynomial 多项式的表达与创建( Expression and Creating of polynomia (1)多项式的表达( expression of polynomial) Maab用行矢量表达多项式系数 Coefficient,各元素按变量的降 幂顺序排列,如多项式为: P(x)=aox"+ajx -+a2x-2.an-1x+an 则其系数矢量 ectorof coefficient)为:P=|a0a1… an-1 an 如将根矢量( lector ofroc)表示为: ar=arar2…arnl 则根矢量与系数矢量之间关系为: (x-ar1)(x-ar2).(x-arn)=aox+ajx -+a2x -2...an-1x+an (2)多项式的创建( polynomial creating) a)系数矢量的直接输入法 利用poly2sym函数直接输入多项式的系数矢量,就可方便的建立 符号形式的多项式。 例:创建多项式x3-4x2+3x+2 poly 2sym([1-43 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+l 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
第四章 MATLAB 的数值计算功能 Chapter 4: Numerical computation of MATLAB 一、多项式(Polynomial)` 1．多项式的表达与创建(Expression and Creating of polynomial) (1) 多项式的表达（expression of polynomial）_ Matlab 用行矢量表达多项式系数(Coefficient)，各元素按变量的降 幂顺序排列，如多项式为： P(x)=a0x n+a1x n-1+a2x n-2…an-1x+an 则其系数矢量(Vector of coefficient)为：P=[a0 a1 … an-1 an] 如将根矢量(Vector of root)表示为： ar=[ ar1 ar2 … arn] 则根矢量与系数矢量之间关系为： （x-ar1）(x-ar2) … (x- arn)= a0x n+a1x n-1+a2x n-2…an-1x+an (2)多项式的创建（polynomial creating） a)系数矢量的直接输入法 利用 poly2sym 函数直接输入多项式的系数矢量，就可方便的建立 符号形式的多项式。 例：创建多项式 x 3 -4x2+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=|638 %根矢量 pa=poly(a) %求系数矢量 ppa=poly2sym(pa) %以符号形式表示原多项式 ezplot(ppa, -50,50D) 1790-144 X^3-17*x^2+90*x-144 注:含复数根的根矢量所创建的多项式要注意: (1)要形成实系数多项式,根矢量中的复数根必须共轭成对; (2)含复数根的根矢量所创建的多项式系数矢量中,可能带有很小的 虚部,此时可采用取实部的命令(rea)把虚部滤掉。 进行多项式的求根运算时,有两种方法,一是直接调用求根函数 roots, poly和 roots互为逆函数。另一种是先把多项式转化为伴随矩阵,然后再 求其特征值,该特征值即是多项式的根。 例3:由给定复数根矢量求多项式系数矢量。 r-0.5-0.3+0.4i-0.3-0.4i; p=poly(r) pr-real(p) ppr=poly2sym(pr)
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)
1.00001.10000.55000.1250 pr 1.00001.10000.55000.1250 ppr- x^3+11/10*x^2+11/20*x+1/8 c)特征多项式输入法 用poly函数可实现由矩阵的特征多项式系数创建多项式。 条件:特征多项式系数矢量的第一个元素必须为一。 例2:求三阶方阵A的特征多项式系数,并转换为多项式形式 a=1638:;756;135 Pa=poly(a) %求矩阵的特征多项式系数矢量 Ppa=poly2sym(pa) 1.0000-160000380000-83.0000 X^3-17x^2+90*x-144 注:n阶方阵的特征多项式系数矢量一定是n+1阶的。 注:(1)要形成实系数多项式,根矢量中的复数根必须共轭成对; (2)含复数根的根矢量所创建的多项式系数矢量中,可能带有很小的 虚部,此时可采用取实部的命令(rea)把虚部滤掉。 进行多项式的求根运算时,有两种方法,一是直接调用求根函数 roots poly和 roots互为逆函数。另一种是先把多项式转化为伴随矩阵,然后再 求其特征值,该特征值即是多项式的根。 例4:将多项式的系数表示形式转换为根表现形式 求x3-6x2-72x-27的根 1-6-72-27 roots(a)
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 3 -6x2 -72x-27 的根 a=[1 -6 -72 -27] r=roots(a)
12.1229 -0.3884 MATLAB约定,多项式系数矢量用行矢量表示,根矢量用列矢量表示。 1.多项式的乘除运算 Multiplication and division of polynomial) 多项式乘法用函数conv(a,b)实现,除法用函数 decoy(a,b)实现。 例1:a(s)=s2+2s+3,b(s)=4s2+5s+6计算a(s)与b(s)的乘积 a=123;b=456; c=conv(a, b) cs=poly2sym(c,’s) 413282718 4*s^4+13*s^3+28*s^2+27*s+18 例2:展开(s2+2s+2)(s+4s+1)(多个多项式相乘) c=conv(1,2,2],conv(1,41,1,1) cs=poly2sym(c,’s %(指定变量为s) 71618 s^4+7s^3+16s^2+18s+8 例2:求多项式S^47*s^3+16^2+18*s+8分别被(s+4),(s+3)除后的结 果 c=1716188 q1r1 decoy(c,1,4)%q—商矢量,r余数矢量 92, r2]=decon(c, [1, 3) cc=conv(q2, 1, 3D) %对除(s+3)结果检验 test=((c-r2)=cc)
r = 12.1229 -5.7345 -0.3884 MATLAB 约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。 >> 1. 多项式的乘除运算(Multiplication and division of polynomial) 多项式乘法用函数 conv(a,b)实现， 除法用函数 deconv(a,b)实现。 例 1：a(s)=s2+2s+3, b(s)=4s2+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： 展开(s2+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)
rI= 0 r2 000 7161818 1.其他常用的多项式运算命令( Other computation command of polynomial) pa= polyval(p,s)按数组运算规则计算给定s时多项式p的值。 pm= polyram(p2s)按矩阵运算规则计算给定s时多项式p的值。 r;p,k= residue(b,a)部分分式展开,b,a分别是分子分母多项式系数 矢量,rp,k分别是留数、极点和直项矢量 p= polyfit(x,y,n)用n阶多项式拟合x,y矢量给定的数据 polder(p) 多项式微分 注:对于多项式b(s)与不重根的n阶多项式a(s)之比,其部分分式 展开为:如=1+5+L+--+(s) a(s) s-p s-p2 S-pm 式中:p1,p2…,pm称为极点( poles),r1,r2y…,rn称为留数( residues), k(s称为直项( direct terms),假如a(s)含有m重根p则相应部分应 写成:++L+ P, (s-P) (s-P, RESIDUE Partial-fraction expansion(residues) R, P, K= RESIDUE(B, A) finds the residues, poles and direct term of a partial raction expansion of the ratio of two polynomials b(SA(s). If there are no multiple
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)之比，其部分分式 展开为： ( ) ( ) ( ) 2 2 1 1 k s s p r L s p r s p r a s b s n n + − + + − + − = 式中：p1,p2,…,pn 称为极点(poles)，r1,r2,…,rn 称为留数(residues)， k(s)称为直项(direct terms)，假如 a(s)含有 m 重根 pj，则相应部分应 写成： m j j m j j j j s p r L s p r s p r ( ) ( ) 1 2 1 − + + − + − + + − 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 R(2) +…+-—+k(s) P(1)s-P(2) 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)+l If Pg=.=P(+m-D)is a pole of multplicity m, then the expansion includes terms of the form R(〔+1) R(〔+m-1) s-P(j(s-P(j)^2(s-P()^ 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:对(3x4+2x3+5x2+4x+6)(x5+3x4+4x3+2x2+7x+2)做部分分式展 开 a=|13427 2l; b=32546 Ir,s, k]=residue(b, a) 1.1274+1.1513i 1.1274-1.1513i 0.0232-007221 0.0232+0.0722i 0.7916
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：对 (3x4+2x3+5x2+4x+6)/(x5+3x4+4x3+2x2+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=12345]; %实验数据 y=5.543.1128290.74984]; p=polyfit(x, y, 3) %做三阶多项式拟合 x2=1:.1:5: 2=polyval(p, x 2); %根据给定值计算多项式结果 pot(x,y,’o’,x2,y2 线性代数( Linear algebra) 解线性方程( Linear equation)就是找出是否存在一个唯一的矩阵x,使 得a,b满足关系: ax=b或xa=b MALAB中x=a\b是方程ax=b的解,x=ba是方程式xa=b的解 通常线性方程多写成ax=b,“”较多用,两者的关系为: (b/a)=(a”b) 系数矩阵a可能是m行n列的,有三种情况: 方阵系统:( Square matrix)m=n可求出精确解(a必须是非奇异 ( nonsingular),即满秩( full rank) 超定系统:( Overdetermined system)m>n可求出最小二乘解 欠定系统:( Underdetermind system)m<n可尝试找出含有最少m个基 解或最小范数解 MATLAB对不同形式的参数矩阵,采用不同的运算法则来处理,它会 自动检测参数矩阵,以区别下面几种形式:
-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) 一． 线性代数(Linear Algebra) 解线性方程(Linear equation)就是找出是否存在一个唯一的矩阵 x，使 得 a,b 满足关系： ax=b 或 xa=b MALAB 中 x=a\b 是方程 ax=b 的解， x=b/a 是方程式 xa=b 的解。 通常线性方程多写成 ax=b，“\”较多用，两者的关系为： (b/a)’=(a’\b’) 系数矩阵 a 可能是 m 行 n 列的，有三种情况： *方阵系统: (Square matrix) m=n 可求出精确解（a 必须是非奇异 (nonsingular)，即满秩(full rank)） *超定系统：(Overdetermind system)m>n 可求出最小二乘解 *欠定系统：(Underdetermind system) m<n 可尝试找出含有最少 m 个基 解或最小范数解 MATLAB 对不同形式的参数矩阵，采用不同的运算法则来处理，它会 自动检测参数矩阵，以区别下面几种形式：
三角矩阵( Triangular matrix) 对称正定矩阵( symmetrical positive determined matrix) 非奇异方阵( Nonsingularmatrix) 超定系统 Overdetermined system) *欠定系统( Underdetermined system) 1.方阵系统:( Square array) 最常见的是系数矩阵为方阵a,常数项b为列矢量,其解x可写成x=ab, x和b大小相同。 例1:求方阵系统的根 a=ll67;5139;1718 b=16134 x=\b 17 b 16 13 3.9763 5.4455 -8.6303 例2:假如a,b为两个大小相同的矩阵,求方阵系统的根 a=459;18195;1413 b=512;31519;7610 x=alb
*三角矩阵(Triangular Matrix) *对称正定矩阵(symmetrical positive determined matrix) *非奇异方阵(Nonsingular matrix) *超定系统(Overdetermind system) *欠定系统(Underdetermind system) 1. 方阵系统：(Square array) 最常见的是系数矩阵为方阵a，常数项b为列矢量, 其解x可写成x=a\b, x 和 b 大小相同。 例 1： 求方阵系统的根。 a=[11 6 7; 5 13 9; 17 1 8] b=[16 13 4]’ x=a\b a = 11 6 7 5 13 9 17 1 8 b = 16 13 4 x = 3.9763 5.4455 -8.6303 例 2：假如 a,b 为两个大小相同的矩阵，求方阵系统的根。 a=[4 5 9; 18 19 5; 1 4 13] b=[1 5 12; 3 15 19; 7 6 10] x=a\b C=a*x
4 1819 12 31519 610 -3.6750-0.733329708 3.725014667-2.1292 -0.32500.06671.1958 1.00005.000012000 3.0000150000190000 7.00006.000010.0000 若方阵a的各个行矢量线性相关 linear correlation),则称方阵a为奇异矩阵。 这时线性方程将有无穷多组解。若方阵是奇异矩阵,则反斜线运算因子将发出 警告信息。 2.超定系统( Overdetermine systen) 实验数据较多,寻求他们的曲线拟合 如在t内测得一组数据y
a = 4 5 9 18 19 5 1 4 13 b = 1 5 12 3 15 19 7 6 10 x = -3.6750 -0.7333 2.9708 3.7250 1.4667 -2.1292 -0.3250 0.0667 1.1958 C = 1.0000 5.0000 12.0000 3.0000 15.0000 19.0000 7.0000 6.0000 10.0000 若方阵a 的各个行矢量线性相关(linear correlation)，则称方阵a 为奇异矩阵。 这时线性方程将有无穷多组解。若方阵是奇异矩阵，则反斜线运算因子将发出 警告信息。 2．超定系统(Overdetermind system) 实验数据较多，寻求他们的曲线拟合。 如在 t 内测得一组数据 y： t y
0.0 0.82 0.3 0.72 0.8 0.63 1.1 0.60 1.6 0.5 2.2 0.50 这些数据显然有衰减指数趋势:y(t)c1+ee 此方程意为y矢量可以由两个矢量逐步逼近而得,一个是单行的常数 矢量,一个是由指数e项构成,两个参数c和c2可用最小二乘法求得,它 们表示实验数据与方程y(t)c1+2e之间距离的最小平方和。 例1:求上述数据的最小二乘解。将数据带入方程式y(t)-c1+c2et中,可 得到含有两个未知数的6个等式,可写成6行2列的矩阵e. t=100.30.81.11.62.2]’; y=0.820.720.630.600.550.50’; e=lones(size(t)) exp(-t)I %求6个y(t)方程的系数矩阵 e=ely %求方程的解 1.00001.0000 1.00000.7408 1.00000.4493 1.00000.3329 1.00000.2019 100000.1108 0.4744 0.3434 带入方程得:y(t)-0.474410.3434e 用此方程可绘制曲线 t=|00.30.81.11.62.2 y=0.820.720.630.600.550.50]’; t1=[0: 0.1: 2.5; yl=lones(size(tD)), exp(-tD)]*c
0.0 0.82 0.3 0.72 0.8 0.63 1.1 0.60 1.6 0.55 2.2 0.50 这些数据显然有衰减指数趋势： y(t)~c1+c2e -t 此方程意为 y 矢量可以由两个矢量逐步逼近而得，一个是单行的常数 矢量，一个是由指数 e -t项构成，两个参数 c1 和 c2 可用最小二乘法求得，它 们表示实验数据与方程 y(t)~c1+c2e -t之间距离的最小平方和。 例 1： 求上述数据的最小二乘解。将数据带入方程式 y(t)~c1+c2e -t 中，可 得到含有两个未知数的 6 个等式，可写成 6 行 2 列的矩阵 e. t=[0 0.3 0.8 1.1 1.6 2.2]’; y=[0.82 0.72 0.63 0.60 0.55 0.50]’; e=[ones(size(t)) exp(-t)] %求 6 个 y(t)方程的系数矩阵 c=e\y % 求方程的解 e = 1.0000 1.0000 1.0000 0.7408 1.0000 0.4493 1.0000 0.3329 1.0000 0.2019 1.0000 0.1108 c = 0.4744 0.3434 带入方程得：y(t)~0.4744+0.3434e-t 用此方程可绘制曲线： t=[0 0.3 0.8 1.1 1.6 2.2]’; y=[0.82 0.72 0.63 0.60 0.55 0.50]’; t1=[0:0.1:2.5]’; y1=[ones(size(t1)),exp(-t1)]*c