[理学]第四章 MATLAB的数值计算功能.doc

第四章 MATLAB 的数值计算功能Chapter 4: Numerical computation of MATLAB 数值计算是MATLAB最基本、最重要的功能，是MATLAB最具代表性的特点。MATLAB在数值计算过程中以数组和矩阵为基础。数组是MATLAB运算中的重要数据组织形式。前面章节对数组、矩阵的特征及其创建与基本运算规则等相关知识已作了较详尽的介绍，本章重点介绍常用的数值计算方法。 一、多项式(Polynomial)多项式在众多学科的计算中具有重要的作用，许多方程和定理都是多项式的形式。MATLAB提供了标准多项式运算的函数，如多项式的求根、求值和微分，还提供了一些用于更高级运算的函数，如曲线拟合和多项式展开等。1多项式的表达与创建(Expression and Creating of polynomial)(1) 多项式的表达（expression of polynomial）_ Matlab用行矢量表达多项式系数(Coefficient)和根，系数矢量中各元素按变量的降幂顺序排列，如多项式为：P(x)=a0xn+a1xn-1+a2xn-2an-1x+an则其系数矢量(Vector of coefficient)为：P=a0 a1 an-1 an如将根矢量(Vector of root)表示为：ar= ar1 ar2 arn则根矢量与系数矢量之间关系为：（x-ar1）(x- ar2) (x- arn)= a0xn+a1xn-1+a2xn-2an-1x+an(2)多项式的创建（polynomial creating）a,系数矢量的直接输入法利用poly2sym函数直接输入多项式的系数矢量，就可方便的建立符号形式的多项式。 例1：创建给定的多项式x3-4x2+3x+2poly2sym(1 -4 3 2)ans =x3-4*x2+3*x+2也可以用poly2str.求一个方阵对应的符号形式的多项式。例 2：a=3 4 6;5 7 1;8 3 9p=poly(a); %求方阵的特征多项式系数矢量poly2str(p,'x') %建立符号形式的多项式ans = x3 - 19 x2 + 40 x + 214POLY 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. POLY2SYM Polynomial coefficient vector to symbolic polynomial.POLY2SYM(C) is a symbolic polynomial in 'x' with coefficients from the vector C.POLY2SYM(C,'V') and POLY2SYM(C,SYM('V') both use the symbolic variable specified by the second argument.POLY2STR Return polynomial as string.S = POLY2STR(P,'s') or S=POLY2STR(P,'z') POLY2SYM Polynomial coefficient vector to symbolic polynomial. POLY2SYM(C,'V') and POLY2SYM(C,SYM('V') both use the symbolic variable specified by the second argument.b) 由根矢量创建多项式多项式行向量可以由命令poly创建，如A为矩阵，则poly（A）将创建A矩阵的特征多项式，如果A为向量 ar1 ar2 arn-1 arn，则创建（x-ar1）(x- ar2) (x-arn-1) (x- arn)生成的多项式的系数矢量。若已知根矢量ar, 通过调用函数 p=poly(ar) 可以产生多项式的系数矢量, 再利用poly2sym函数就可方便的建立符号形式的多项式。例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 -144ppa =x3-17*x2+90*x-144说明：（1）根矢量元素为n ，则多项式系数矢量元素为n+1； （2）函数poly2sym(pa) 把多项式系数矢量表达成符号形式的多项式，缺省情况下自变量符号为x，可以指定自变量。（3）使用简单绘图函数ezplot可以直接绘制符号形式多项式的曲线。例 2： 由给定复数根矢量求多项式系数矢量。r=-0.5 -0.3+0.4i -0.3-0.4i;p=poly(r) %求系数矢量pr=real(p) %提取系数矢量实部 ppr=poly2sym(pr) %以符号形式表示原多项式p = 1.0000 1.1000 0.5500 0.1250pr = 1.0000 1.1000 0.5500 0.1250ppr =x3+11/10*x2+11/20*x+1/8说明：含复数根的根矢量所创建的多项式要注意： （1）要形成实系数多项式，根矢量中的复数根必须共轭成对； （2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令real把虚部滤掉。c) 多项式的求根如果需要进行系数表示形式的多项式的求根运算，有两种方法可以实现，一是直接调用求根函数roots，poly和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。d) 特征多项式输入法用poly函数可实现由矩阵的特征多项式系数创建多项式。条件：特征多项式系数矢量的第一个元素必须为1。例 1： 求三阶方阵A的特征多项式系数，并转换为多项式形式。a=6 3 8;7 5 6; 1 3 5Pa=poly(a) %求矩阵的特征多项式系数矢量Ppa=poly2sym(pa) %建立特征多项式Pa = 1.0000 -16.0000 38.0000 -83.0000Ppa =x3-17*x2+90*x-144注：n 阶方阵的特征多项式系数矢量一定是n +1阶的。例 2： 将多项式的系数表示形式转换为根表现形式（求根）。求 x3-6x2-72x-27的根a=1 -6 -72 -27r=roots(a)r = 12.1229 -5.7345 -0.3884MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。2. 多项式的运算(Computation of polynomial)a, 多项式的加减法运算(addition and subtraction of polynomial)MATLAB没有提供多项式直接的加减法函数，若两个多项式向量大小相同，标准的数组加法有效，当两个多项式阶次不同时，低阶多项式必须用首零填补使其与高阶多项式有相同的阶次。例1：a=3 2 5 8;b=3 5 1 9; c=a+b %同价求和 d=2 4 1 -4 5; e=0 c+d %不同价求和c = 6 7 6 17e = 2 10 8 2 22减法的操作与加法类似，相当于将后项变量取反。b,多项式的乘除运算(Multiplication and division of polynomial)多项式的乘除法对应着向量的卷积与解卷积，多项式乘法（卷积）用函数conv(a,b)实现， 除法（解卷积）用函数deconv(a,b)实现。长度为m的向量a和长度为n的向量b的卷积定义为：C(k)=C向量的长度为：m+n-1解卷积是卷积的逆运算，向量a对向量c进行解卷积将得到商向量 q和余量r，并且满足：例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) %建立指定变量为s的符号形式多项式c = 4 13 28 27 18cs =4*s4+13*s3+28*s2+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 8cs =s4+7*s3+16*s2+18*s+8例3：求多项式s4+7*s3+16*s2+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)q1 = 1 3 4 2r1 = 0 0 0 0 0q2 = 1 4 4 6r2 = 0 0 0 0 -10cc = 1 7 16 18 18test = 1 1 1 1 13. 其他常用的多项式运算命令(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分别是分子（numerator）分母（denominator）多项式系数矢量，r,p,k分别是留数、极点和直项矢量p=polyfit(x,y,n) 用n阶多项式拟合x，y矢量给定的数据。polyder(p) 多项式微分。注： 对于多项式b(s)与不重根的n阶多项式a(s)之比，其部分分式展开为： 式中：p1,p2,pn称为极点(poles)，r1,r2,rn 称为留数(residues)，k(s)称为直项(direct terms)，假如a(s)含有m重根pj，则相应部分应写成：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,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 residuesare 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)mB,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.例1：对 (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将是空矩阵，表示无此项)例2：对一组实验数据进行多项式最小二乘拟合(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=bMALAB中左除x=ab 是方程式 ax=b 的解，右除b/a是方程式xa=b的解。通常线性方程多写成ax=b，“”较多用，左、右除两者的关系为：(b/a)=(ab)系数矩阵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)*非奇异方阵(Nonsingular matrix)*超定系统(Overdetermind system)*欠定系统(Underdetermind system)1. 方阵系统：(Square array) m=n最常见的方程式为ax=b形式，系数矩阵a为方阵，常数项b为列矢量, 其解x可写成x=ab, x和b大小相同。例1： 用左除求方阵系统的根。a=11 6 7; 5 13 9; 17 1 8b=16 13 4x=aba = 11 6 7 5 13 9 17 1 8b = 16 13 4x = 3.9763 5.4455 -8.6303例2：假如a,b为两个大小相同的矩阵，可用左除求方阵系统的根。a=4 5 9; 18 19 5; 1 4 13b=1 5 12; 3 15 19; 7 6 10x=abC=a*xa = 4 5 9 18 19 5 1 4 13b = 1 5 12 3 15 19 7 6 10x = -3.6750 -0.7333 2.9708 3.7250 1.4667 -2.1292 -0.3250 0.0667 1.1958C = 1.0000 5.0000 12.0000 3.0000 15.0000 19.00007.0000 6.0000 10.0000若方阵a的各个行矢量线性相关(linear correlation)，则称方阵a为奇异矩阵。这时线性方程将有无穷多组解。若方阵是奇异矩阵，则反斜线运算因子将发出警告信息。2超定系统(Overdetermind system) mn 实验数据较多，寻求他们的曲线拟合。 如在t内测得一组数据y：t y0.0 0.820.3 0.720.8 0.631.1 0.601.6 0.552.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=ey % 求方程的解e = 1.0000 1.0000 1.0000 0.7408 1.0000 0.4493 1.0000 0.3329 1.0000 0.2019 1.0000 0.1108c = 0.47440.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)*cplot(t1,y1,b,t,y,ro)如果一个矩阵的行矢量是线性相关的，则它的最小二乘解并不唯一，因此，ab运算将给出警告，并产生含有最少元素的基解。3 .欠定系统: (Underdetermind system) mn欠定系统为线性相关系统，其解都不唯一，MATLAB会计算一组构成通解的基解，而方程的特解则用QR分解法决定。两种解法：最少元素解ab，最小范数解pinv(a)*b.例1： 用两种方法求解欠定系统。对a和矢量b分别用ab和pinv(a)*b求解:a=1 1 1; 1 1 -1b=10 6p=abq=pinv(a)*ba = 1 1 1 1 1 -1b = 10 6p = 8.0000 0 2.0000q = 4.0000 4.0000 2.0000三 逆矩阵及行列式(Revers and determinant of matrix)1 方阵的逆和行列式(Revers and determinant of square matrix)若a是方阵，且为非奇异阵，则方程ax=I和 xa=I有相同的解X。X称为a的逆矩阵，记做a-1，在MATLAB中 用inv 函数来计算矩阵的逆。计算方阵的行列式则用det函数。DET Determinant.DET(X) is the determinant of the square matrix X. Use COND instead of DET to test for matrix singularity.INV Matrix inverse.INV(X) is the inverse of the square matrix X. A warning message is printed if X is badly scaled or nearly singular.例1：计算方阵的行列式和逆矩阵。a=3 -3 1;-3 5 -2;1 -2 1;b=14 13 5; 5 1 12;6 14 5;d1=det(a) %求方阵的行列式x1=inv(a) %求逆d2=det(b)x2=inv(b)d1 = 1x1 = 1.0000 1.0000 1.0000 1.0000 2.0000 3.0000 1.0000 3.0000 6.0000d2 = -1351x2 = 0.1207 -0.0037 -0.1118 -0.0348 -0.0296 0.1058 -0.0474 0.0873 0.03772 广义逆矩阵（伪逆）(Generalized inverse matrix)一般非方阵无逆矩阵和行列式，方程ax=I 和xa=I至少有一个无解，这种矩阵可以求得特殊的逆矩阵，称为广义逆矩阵(generalized inverse matrix)（或伪逆 pseudoinverse）。矩阵amn存在广义逆矩阵xnm，使得 ax=Imn， MATLAB用pinv函数来计算广义逆矩阵。例1：计算广义逆矩阵。a=8 14; 1 3; 9 6x=pinv(a)b=x*ac=a*xd=c*a %d=a*x*a=ae=x*c %e=x*a*x=xa = 8 14 1 3 9 6x = -0.0661 -0.0402 0.1743 0.1045 0.0406 -0.0974b = 1.0000 -0.0000 -0.0000 1.0000c = 0.9334 0.2472 0.0317 0.2472 0.0817 -0.1177 0.0317 -0.1177 0.9849d = 8.0000 14.0000 1.0000 3.0000 9.0000 6.0000e = -0.0661 -0.0402 0.17430.1045 0.0406 -0.0974PINV Pseudoinverse.X = PINV(A) produces a matrix X of the same dimensions as A' so that A*X*A = A, X*A*X = X and A*X and X*A are Hermitian. The computation is based on SVD(A) and any singular values less than a tolerance are treated as zero. The default tolerance is MAX(SIZE(A) * NORM(A) * EPS.PINV(A,TOL) uses the tolerance TOL instead of the default.四 矩阵分解(Matrix decomposition) 通过矩阵分解的方法求解大型方程组非常有效的，这种方法可以使运算速度加快，节省磁盘空间，节省内存。MATLAB求解线性方程的过程基于三种分解法则：（）Cholesky分解，针对对称正定矩阵；（）LU分解，高斯消元法， 针对一般矩阵； （）QR分解，正交化， 针对一般长方形矩阵（行数列数）这三种分解运算分别由chol, lu和 qr三个函数来分解.1 Cholesky分解(Cholesky Decomposition)仅适用于对称和上三角矩阵，如果A为对称正定矩阵，则Cholesky分解可将矩阵A分解为上三角矩阵和其转置的乘积， 即：A=R*R, R为上三角阵。方程A*X=b变成R*R*X=b, 所以有 X=R(Rb)例1：cholesky分解。a=pascal(6)b=chol(a)a = 1 1 1 1 1 1 1 2 3 4 5 6 1 3 6 10 15 21 1 4 10 20 35 56 1 5 15 35 70 126 1 6 21 56 126 252b = 1 1 1 1 1 1 0 1 2 3 4 5 0 0 1 3 6 10 0 0 0 1 4 10 0 0 0 0 1 50 0 0 0 0 1CHOL Cholesky factorization.CHOL(X) uses only the diagonal and upper triangle of X. The lower triangular is assumed to be the (complex conjugate) transpose of the upper. If X is positive definite, then R = CHOL(X) produces an upper triangular R so that R'*R = X. If X is not positive definite, an error message is printed.R,p = CHOL(X), with two output arguments, never produces anerror message. If X is positive definite, then p is 0 and R is the same as above. But if X is not positive definite, then p is a positive integer. When X is full, R is an upper triangular matrix of order q = p-1so that R'*R = X(1:q,1:q). When X is sparse, R is an upper triangular matrix of size q-by-n so that the L-shaped region of the first q rows and first q columns of R'*R agree with those of X.例2：对下列矩阵进行cholesky分解，并求解方程组的解。16x1+4x2+8x3=284x1+5x2-4x3=58x1-4x2+22x3=26A=16 4 8;4 5 -4;8 -4 22;b=28 5 26'R=chol(A) %Cholesky分解X=R(R'b) %根据R矩阵求解方程组的解R = 4 1 2 0 2 -3 0 0 3X = 1 1 12 LU分解(LU factorization).用lu函数完成LU分解，将矩阵分解为上、下两个三角阵，其调用格式为：l,u=lu(a) l代表下三角阵，u代表上三角阵。例1：LU分解。a=47 24 22; 11 44 0;30 38 41l,u=lu(a)a = 47 24 22 11 44 0 30 38 41l = 1.0000 0 0 0.2340 1.0000 0 0.6383 0.5909 1.0000u = 47.0000 24.0000 22.0000 0 38.3830 -5.1489 0 0 30.0000LU LU factorization.L,U = LU(X) stores an upper triangular matrix in U and a "psychologically lower triangular matrix" (i.e. a product of lower triangular and permutation matrices) in L, so that X = L*U. X can be rectangular. L,U,P = LU(X) returns unit lower triangular matrix L, upper triangular matrix U, and permutation matrix P so that P*X = L*U.3 QR分解(Orthogonal-triangular decomposition).函数调用格式：q,r=qr(a), q代表正规正交矩阵，r代表三角形矩阵。原始阵a不必一定是方阵。如果矩阵a是m×n阶的，则矩阵q是m×m阶的，矩阵r是m×n阶的。例1：QR分解.A=22 46 20 20; 30 36 46 44;39 8 45 2;q,r=qr(A)q = -0.408