1、第四章 MATLAB 的数值计算功能Chapter 4: Numerical computation of MATLAB数值计算是 MATLAB 最基本、最重要的功能，是 MATLAB 最具代表性的特点。MATLAB 在数值计算过程中以数组和矩阵为基础。数组是MATLAB 运算中的重要数据组织形式。前面章节对数组、矩阵的特征及其创建与基本运算规则等相关知识已作了较详尽的介绍，本章重点介绍常用的数值计算方法。一、多项式(Polynomial)多项式在众多学科的计算中具有重要的作用，许多方程和定理都是多项式的形式。MATLAB 提供了标准多项式运算的函数，如多项式的求根、求值和微分，还提供了一些用
2、于更高级运算的函数，如曲线拟合和多项式展开等。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-ar
3、1）(x- ar 2) (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+2POLY Convert roots to polynomial.POLY(A), when A is an N by N matrix, is a row vector withN+1 elements which a
4、re the coefficients of the characteristic polynomial, DET(lambda*EYE(SIZE(A) - A) .POLY(V), when V is a vector, is a vector whose elements arethe coefficients of the polynomial whose roots are theelements of V . For vectors, ROOTS and POLY are inversefunctions of each other, up to ordering, scaling,
5、 androundoff error.b) 由根矢量创建多项式已知根矢量 ar, 通过调用函数 p=poly(ar)产生多项式的系数矢量, 再利用 poly2sym 函数就可方便的建立符号形式的多项式。例 2：由根矢量创建多项式。将多项式(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）函数 p
6、oly2sym(pa) 把多项式系数矢量表达成符号形式的多项式，缺省情况下自变量符号为 x，可以指定自变量。（3）使用简单绘图函数 ezplot 可以直接绘制符号形式多项式的曲线。例 3： 由给定复数根矢量求多项式系数矢量。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）要形成实系数多项式，根矢
7、量中的复数根必须共轭成对；（2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real) 把虚部滤掉。如果需要进行系数表示形式的多项式的求根运算，有两种方法可以实现，一是直接调用求根函数 roots，poly 和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。c) 特征多项式输入法用 poly 函数还可实现由矩阵的特征多项式系数创建多项式。条件：特征多项式系数矢量的第一个元素必须为 1。例 2： 求三阶方阵 A 的特征多项式系数，并转换为多项式形式。a=6 3 8;7 5 6; 1 3 5Pa=pol
8、y(a) %求矩阵的特征多项式系数矢量Ppa=poly2sym(pa)Pa =1.0000 -16.0000 38.0000 -83.0000Ppa =x3-17*x2+90*x-144注：n 阶方阵的特征多项式系数矢量一定是 n +1 阶的。例 4： 将多项式的系数表示形式转换为根表现形式。求 x3-6x2-72x-27 的根a=1 -6 -72 -27r=roots(a)r =12.1229-5.7345-0.3884MATLAB 约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。2. 多项式的乘除运算(Multiplication and division of polynomial
9、)向量的卷积与解卷积对应着多项式的乘除法，多项式乘法（卷积）用函数 conv(a,b)实现， 除法（解卷积）用函数 deconv(a,b)实现。长度为 m 的向量 a 和长度为 n 的向量 b 的卷积定义为：C(k)=1()(+1)C 向量的长度为：m+n-1解卷积是卷积的逆运算，向量 a 对向量 c 进行解卷积将得到商向量 q 和余量 r，并且满足：（ ） -r(k)=1()(+1)例 1：a(s)=s 2+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) %建
10、立指定变量为 s 的符号形式多项式c =4 13 28 27 18cs =4*s4+13*s3+28*s2+27*s+18例 2： 展开(s 2+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,
11、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 分别是分子
12、（numerator）分母（denominator）多项式系数矢量，r,p,k 分别是留数、极点和直项矢量p=polyfit(x,y,n) 用 n 阶多项式拟合 x，y 矢量给定的数据。polyder(p) 多项式微分。注： 对于多项式 b(s)与不重根的 n 阶多项式 a(s)之比，其部分分式展开为： )()(21 skpsrLsrpsabn式中：p 1,p2,pn 称为极点 (poles)，r 1,r2,rn 称为留数(residues)，k(s)称为直项(direct terms)，假如 a(s)含有 m 重根 pj， 则相应部分应写成： mjjj psrLpsr)()(121RESID
13、UE 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
14、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 te
15、rm coefficient vector is empty if length(B) n 可求出最小二乘解*欠定系统：(Underdetermind system) m n, then only the first n columns of Q are computed.4. 特征值与特征矢量(Eigenvalues and eigenvectors).MATLAB 中使用函数 eig 计算特征值和 特征矢量，有两种调用方法：*e=eig(a), 其中 e 是包含特征值的矢量；*v,d=eig(a), 其中 v 是一个与 a 相同的 nn 阶矩阵，它的每一列是矩阵 a 的一个特征值所对应的特
16、征矢量，d 为对角阵，其对角元素即为矩阵a 的特征值。例：计算特征值和特征矢量。a=34 25 15; 18 35 9; 41 21 9e=eig(a)v,d=eig(a)a =34 25 1518 35 941 21 9e =68.506615.5122-6.0187v =-0.6227 -0.4409 -0.3105-0.4969 0.6786 -0.0717-0.6044 -0.5875 0.9479d =68.5066 0 00 15.5122 00 0 -6.0187EIG Eigenvalues and eigenvectors.E = EIG(X) is a vector con
17、taining the eigenvalues of a square matrix X.V,D = EIG(X) produces a diagonal matrix D of eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that X*V = V*D.V,D = EIG(X,nobalance) performs the computation with balancingdisabled, which sometimes gives more accurate res
18、ults for certainproblems with unusual scaling. If X is symmetric, EIG(X,nobalance)is ignored since X is already balanced.5. 奇异值分解.( Singular value decomposition).如存在两个矢量 u,v 及一常数 c,使得矩阵 A 满足：Av=cu, Au=cv称 c 为奇异值，称 u,v 为奇异矢量。将奇异值写成对角方阵，而相对应的奇异矢量作为列矢量则可写成两个正交矩阵 U，V， 使得： AV=U ， AU=V 因为 U，V 正交，所以可得奇异值表达
19、式：A=UV 。一个 m 行 n 列的矩阵 A 经奇异值分解，可求得 m 行 m 列的 U, m 行 n列的矩阵和 n 行 n 列的矩阵 V.。奇异值分解用 svd 函数实现，调用格式为；u,s,v=svd(x) SVD Singular value decomposition.U,S,V = SVD(X) produces a diagonal matrix S, of the same dimension as X and with nonnegative diagonal elements in decreasing order, and unitary matrices U and V
20、 so that X = U*S*V.S = SVD(X) returns a vector containing the singular values.U,S,V = SVD(X,0) produces the “economy size“ decomposition. If X is m-by-n with m n, then only the first n columns of U are computed and S is n-by-n.例： 奇异值分解。a=8 5; 7 3;4 6;u,s,v=svd(a) % s 为奇异值对角方阵u =-0.6841 -0.1826 -0.70
21、61-0.5407 -0.5228 0.6591-0.4895 0.8327 0.2589s =13.7649 00 3.08650 0v =-0.8148 -0.5797-0.5797 0.8148五 数据分析(Data Analysis)MATLAB 对数据分析有两条约定：（1） 若输入量 X 是矢量，则不论是行矢量还是列矢量，运算是对整个矢量进行的；（2）若输入量 X 是数组， （或称矩阵） ，则命令运算是按列进行的。即默认每个列是有一个变量的不同“观察“所得的数据组成。 1. 基本统计命令 （表 4-1）例： 做各种基本统计运算。A=5 -10 -6 0;2 6 3 -3;-9 5 -
22、10 11;-22 17 0 -19;-1 6 -4 4Amax=max(A) %找 A 各列的最大元素Amin=min(A) %找 A 各列的最小元素Amed=median(A) %列数为奇数使找 A 各列中数值居中的中位元素，列数为偶数使取两个中位元素的平均值Amean=mean(A) %找 A 各列的平均值Astd=std(A) %求 A 各列的标准差Aprod=prod(A) %求 A 各列元素的积Asum=sum(A) %求 A 各列元素的和S=cumsum(A) %求 A 各列元素的累积和P=cumprod(A) %求 A 各列元素的累积 j 积I=sort(A) %使 A 的各列
23、元素按递增排列A =5 -10 -6 02 6 3 -3-9 5 -10 11-22 17 0 -19-1 6 -4 4Amax =5 17 3 11Amin =-22 -10 -10 -19Amed =-1 6 -4 0Amean =-5.0000 4.8000 -3.4000 -1.4000Astd =10.8397 9.6281 5.0794 11.1490Aprod =-1980 -30600 0 0Asum =-25 24 -17 -7S =5 -10 -6 07 -4 -3 -3-2 1 -13 8-24 18 -13 -11-25 24 -17 -7P =5 -10 -6 010
24、 -60 -18 0-90 -300 180 01980 -5100 0 0-1980 -30600 0 0I =-22 -10 -10 -19-9 5 -6 -3-1 6 -4 02 6 0 45 17 3 11求矩阵元素的最大值、最小值可用：Amax=max(max(A) 或 Amax=max(A(:)， Amin=min(min(A) 或 Amin=min(A(:) 2协方差阵和相关阵(Covariance matrix and Correlation coefficients).向量的协方差体现了向量中各元素的分散程度，矩阵的协方差体现了矩阵中各列元素的相关性。求协方差阵和相关阵的命令
25、：C=cov(x) 求协方差阵C=cov(x,y) 求两随机变量的协方差P=corrcoef(x) 关系数函数，是向量间相关性的归一化表示，求相关阵P=corrcoef(x,y) 求两随机变量的（ 22）相关系数,是向量间相关性的归一化表示。这些函数应用于 m 行 n 列矩阵时结果均为 n 行 n 列的对称矩阵，其第p 行、第 q 列的元素表示矩阵的第 p 行、第 q 列的协方差或相关系数。例 1： 计算协方差和相关阵。x=rand(10,3); y=rand(10,3);cx=cov(x) %求协方差阵cy=cov(y)cxy=cov(x,y) %求两随机变量的协方差px=corrcoef(
26、x) %求相关阵pxy=corrcoef(x,y) %求两随机变量的（22）相关系数cx =0.0483 -0.0066 0.0146-0.0066 0.0283 0.01540.0146 0.0154 0.0978cy =0.1177 0.0073 -0.01270.0073 0.0239 -0.0230-0.0127 -0.0230 0.0772cxy =0.0550 0.00230.0023 0.0697px =1.0000 -0.1783 0.2118-0.1783 1.0000 0.29340.2118 0.2934 1.0000pxy =1.0000 0.03720.0372 1.
27、0000COV Covariance matrix.COV(X), if X is a vector, returns the variance. For matrices, where each row is an observation, and each column a variable, COV(X) is the covariance matrix. DIAG(COV(X) is a vector of variances for each column, and SQRT(DIAG(COV(X) is a vector of standard deviations. COV(X,
28、Y), where X and Y are vectors of equal length, is equivalent to COV(X(:) Y(:). COV(X) or COV(X,Y) normalizes by (N-1) where N is the number of observations. This makes COV(X) the best unbiased estimate of thecovariance matrix if the observations are from a normal distribution.CORRCOEF Correlation co
29、efficients.R=CORRCOEF(X) calculates a matrix R of correlation coefficients for an array X, in which each row is an observation and each column is a variable.R=CORRCOEF(X,Y), where X and Y are column vectors, is the same asR=CORRCOEF(X Y). If C is the covariance matrix, C = COV(X), then CORRCOEF(X) i
30、s the matrix whose (i,j)th element isC(i,j)/SQRT(C(i,i)*C(j,j).2、求极值：在数学上，可以通过确定函数导数为零的点，求出极值点。因为 f(x)的极小值问题等价于-f(x)的极大值问题, 所以 MATLAB 的 Toolbox 函数中只有处理极小值命令的函数。 MATLAB 提供了 3 个从数值上寻找函数极值点的函数。x,fual,exitflag,output=fminbnd(fun,x1,x2,options): 求一元函数在区间(x1,x2)中极小值的指令的最完整格式。x,fual,exitflag,output=fminsea
31、rch(fun,x0,options): 单纯形法求多元函数极值点指令的最完整格式。x,fual,exitflag,output，grad，hessian=fminunc(fun,x0,options): 拟牛顿法求多元函数极值点指令的最完整格式。其中 x1,x2 分别表示被研究区间的左、右边界。 fminsearch 和fminunc 中的 x0 是一个向量，表示极值点位置的初步推测。输入变量options 的默认值可以用 options=optimset(FunFun_Name)获得。例，求解函数的极值%求最小值fn=2*exp(-x)*sin(x); %定义 fn 函数xmin=fmin
32、bnd(fn,2,5) %在区间（2， 5） 寻找最小值emin=5*pi/4-xmin %求最小值的误差x=xmin; %令 x 为最小值点ymin=eval(fn) %计算最小值点的函数值%求最大值fx=-(2*exp(-x)+sin(x); % 为求最大值，函数值取负xmax=fminbnd(fx,0,3) %在区间（0，3）寻找最大值emax=pi/4-xmax %求最大值的误差ymax=eval(fx) %求最大值点的函数值单纯形法求多元函数极值点。fx=2*exp(-x)*sin(x);x=fminsearch(fx,1) %用单纯法求解，从点 1 开始搜索3、数值积分：数值积分函
33、数quad(fname,a,b,tol,trace,p1,p2,) 低阶法-自适应递推辛普生法quadl(fname,a,b,tol,trace,p1,p2,) 高阶法- 自适应递推牛顿-柯西法trapz（x,y） 计算变量为 x，函数 y 的积分fname 是被积函数表达式字符串或函数文件名，a,b 分别为积分的上、下限，tol 用来控制积分精度，默认为 tol=0.001,trace 取 1，表示用图形表示积分结果，取 0 则无图形，p1,p2是向函数传递的参数，可取默认值。例：用 trapz 在区间-1x3 上计算 y=humps(x)下面的面积x=-1:0.12:3;y=humps(x
34、);% Y = HUMPS(X) is a function with strong maxima near % x = .3 and x = .9.plot(x,y) area=trapz(x,y)%以更好的精度计算x=-1:0.02:3;y=humps(x);area=trapz(x,y)%用 quad 和 quadl 在区间- 1x3 上计算 y=humps(x)下面的面积format longarea=quad(humps,-1,3)area=quadl(humps,-1,3)4. 数值微分与梯度 (Difference and approximate derivative ,grad
35、ient).数值微分函数 diff(x)例 1：按列求微分。x=1,10,20;2,12,23;3,14,26;3,16,29d=diff(x)figure(1)plot(x)figure(2)plot(d) %求一阶微分x =1 10 202 12 233 14 263 16 29d =1 2 31 2 30 2 3数值微分比较困难，函数的微小变化都容易产生相邻斜率的大的变化。因而应尽可能避免数值微分。特别是对实验数据进行微分的时候，最好用最小二次曲线拟合数据，然后对所得的多项式进行微分。 x=0:0.1:1;y=-0.443 1.783 3.231 6.21 7.08 7.73 7.98
36、8.99 9.23 9.93 10.99 ;n=2;p=polyfit(x,y,n) %多项式拟合xi=linspace(0,1,100);z=polyval(p,xi);figure(1)hold onplot(x,y,o,x,y,xi,z,:)xlabel(x)ylabel(y)title(二次曲线拟合 )pd=polyder(p) %多项式微分w=polyval(pd,xi);plot(xi,w,*)xlabel(x)ylabel(y)title(曲线微分 )%函数 diff 是对一些描述某函数的数据进行粗略计算的微分函数，用于计算数组中元素间的差分。f=f(x) 的微分近似为 dy/d
37、x=f(x+h)-f(x)/(x+h)-xdy=diff(y)./diff(x); xd=x(1:max(size(x)-1);figure(2)plot(xd,dy)title(用 diff 进行近似差分)ylabel(dy/dx 轴)xlabel(x 轴)结果可见用有限差分近似微分导致很差的结果。例 2： 对于（u=x 2+y2 和 2=4）求 5 点差分。x,y=meshgrid(-4:4,-3:3)u=x.2+y.2v4=4*del2(u) %求 mn 阶矩阵 U 的五点差分矩阵u =25 18 13 10 9 10 13 18 2520 13 8 5 4 5 8 13 2017 10
38、 5 2 1 2 5 10 1716 9 4 1 0 1 4 9 1617 10 5 2 1 2 5 10 1720 13 8 5 4 5 8 13 2025 18 13 10 9 10 13 18 25v4 =4 4 4 4 4 4 4 4 44 4 4 4 4 4 4 4 44 4 4 4 4 4 4 4 44 4 4 4 4 4 4 4 44 4 4 4 4 4 4 4 44 4 4 4 4 4 4 4 44 4 4 4 4 4 4 4 4MESHGRID X and Y arrays for 3-D plots.X,Y = MESHGRID(x,y) transforms the do
39、main specified by vectorsx and y into arrays X and Y that can be used for the evaluationof functions of two variables and 3-D surface plots.The rows of the output array X are copies of the vector x andthe columns of the output array Y are copies of the vector y.X,Y = MESHGRID(x) is an abbreviation f
40、or X,Y = MESHGRID(x,x).X,Y,Z = MESHGRID(x,y,z) produces 3-D arrays that can be used toevaluate functions of three variables and 3-D volumetric plots.DEL2 Discrete Laplacian.L = DEL2(U) when U is a matrix, is an discrete approximation of 0.25*del2 u = (d2u/dx2 + d2/dy2)/4. The matrix L is the same si
41、ze as U with each element equal to the difference between an element of U and the average of its four neighbors.L = DEL2(U) when U is an N-D array, returns an approximation of(del2 u)/2/n where n is ndims(u).L = DEL2(U,H), where H is a scalar, uses H as the spacing betweenpoints in each direction (H
42、=1 by default).L = DEL2(U,HX,HY) when U is 2-D, uses the spacing specified by HXand HY. If HX is a scalar, it gives the spacing between points inthe x-direction. If HX is a vector, it must be of length SIZE(U,2)and specifies the x-coordinates of the points. Similarly, if HYis a scalar, it gives the
43、spacing between points in they-direction. If HY is a vector, it must be of length SIZE(U,1) andspecifies the y-coordinates of the points.L = DEL2(U,HX,HY,HZ,.) when U is N-D, uses the spacing given byHX, HY, HZ, etc.例 3：产生一个二元函数偏导数和梯度。x=-2:0.2:2;y=-2:0.2:2;xx,yy=meshgrid(x,y);z=xx.*exp(-xx.2-yy.2);G
44、x,Gy=gradient(z,0.2,0.2); %Gx,Gy 分别是二元函数的偏导contour(x,y,z,k),hold on,quiver(xx,yy,Gx,Gy,r),hold offX,Y = MESHGRID(x,y) transforms the domain specified by vectorsx and y into arrays X and Y that can be used for the evaluation of functions of two variables and 3-D surface plots. The rows of the output
45、array X are copies of the vector x and the columns of the output array Y are copies of the vector y.DIFF Difference and approximate derivative.DIFF(X), for a vector X, is X(2)-X(1) X(3)-X(2) . X(n)-X(n-1).DIFF(X), for a matrix X, is the matrix of row differences,X(2:n,:) - X(1:n-1,:).DIFF(X), for an
46、 N-D array X, is the difference along the firstnon-singleton dimension of X.DIFF(X,N) is the N-th order difference along the first non-singleton dimension (denote it by DIM). If N = size(X,DIM), DIFF takes successive differences along the next non-singleton dimension.（FX,FY = GRADIENT(F,HX,HY), when F is 2-D, uses the spacingspecified by HX and HY. HX and HY can either be scalars to specifythe spacing between coordinates or vectors to specify the coordinates of the points. If HX and HY are
