第5章-MATLAB数值计算教学课件

5.1 特殊矩阵特殊矩阵5.2 矩阵分析矩阵分析5.3 矩阵分解与线性方程组求解矩阵分解与线性方程组求解5.4 数据处理与多项式计算数据处理与多项式计算5.5 傅立叶分析傅立叶分析5.6 数值微积分数值微积分5.7 常微分方程的数值求解常微分方程的数值求解5.8 非线性方程的数值求解非线性方程的数值求解5.9 稀疏矩阵稀疏矩阵5.1 特殊矩阵特殊矩阵5.1.1对角阵与三角阵对角阵与三角阵1.矩阵的对角元素矩阵的对角元素(1)提取矩阵的对角线元素提取矩阵的对角线元素 设设A为为mn矩阵，矩阵，diag(A)函数用于提取矩阵函数用于提取矩阵A主对角线元主对角线元素产生一个具有素产生一个具有min(m,n)个元素的列向量。个元素的列向量。diag(A)函数还有更进一步的形式函数还有更进一步的形式diag(A,k)，其功能是提，其功能是提取第取第k条对角线的元素。条对角线的元素。目录目录A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19A=17 0 1 0 15 23 5 7 14 16 4 0 13 0 22 10 12 19 21 3 11 18 25 2 19 diag(A)ans=17 5 13 21 19 diag(A,3)ans=0 16(2)构造对角矩阵构造对角矩阵 设设V为具有为具有m个元素的向量，个元素的向量，diag(V)将产生一将产生一个个mm对角矩阵，其主对角线元素即为向量对角矩阵，其主对角线元素即为向量V的的元素。元素。diag(V)函数也有更进一步的形式函数也有更进一步的形式diag(V,k)，其功能是产生一个其功能是产生一个nn(n=m+)对角阵，其第对角阵，其第k条条对角线的元素即为向量对角线的元素即为向量V的元素。的元素。V=1 2 3 4 5;diag(V)ans=1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 4 0 0 0 0 0 5diag(V,2)ans=0 0 1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 例例5.1 先建立先建立55矩阵矩阵A，然后将，然后将A的第的第1行元素乘行元素乘以以1，第，第2行乘以行乘以2，第，第5行乘以行乘以5。目录目录ans=17 0 1 0 15 46 10 14 28 32 12 0 39 0 66 40 48 76 84 12 55 90 125 10 95A=17 0 1 0 15 23 5 7 14 16 4 0 13 0 22 10 12 19 21 3 11 18 25 2 19命令如下：命令如下：A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19;D=diag(1,2,3,4,5);D*A 2.矩阵的三角阵矩阵的三角阵 (1)下三角矩阵下三角矩阵 求矩阵求矩阵A的下三角阵的的下三角阵的MATLAB函数是函数是tril(A)tril(A)函数也有更进一步的一种形式函数也有更进一步的一种形式tril(A,k)，其功能是求矩阵其功能是求矩阵A的第的第k条对角线以下的元素。条对角线以下的元素。(2)上三角矩阵上三角矩阵 在在MATLAB中，提取矩阵中，提取矩阵A的上三角矩阵的函的上三角矩阵的函数是数是triu(A)和和triu(A,k)，其用法与提取下三角，其用法与提取下三角矩阵的函数矩阵的函数tril(A)和和tril(A,k)完全相同。完全相同。目录目录 tril(A)ans=17 0 0 0 0 23 5 0 0 0 4 0 13 0 0 10 12 19 21 0 11 18 25 2 19 triu(A)ans=17 0 1 0 15 0 5 7 14 16 0 0 13 0 22 0 0 0 21 3 0 0 0 0 19 A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19A=17 0 1 0 15 23 5 7 14 16 4 0 13 0 22 10 12 19 21 3 11 18 25 2 19 tril(A,1)ans=17 0 0 0 0 23 5 7 0 0 4 0 13 0 0 10 12 19 21 3 11 18 25 2 19 5.1.2 特殊矩阵的生成特殊矩阵的生成 1.魔方矩阵魔方矩阵魔方矩阵是魔方矩阵是n*n元素所构成的方阵，其每个元素元素所构成的方阵，其每个元素由不同的由不同的1n2的整数所组成，它的每行、每列的整数所组成，它的每行、每列以及对角线元素之和均相等，并等于以及对角线元素之和均相等，并等于n(1+n2)/2.函数格式为函数格式为magic(n)例例5.2 将将101125等等25个数填入一个个数填入一个5行行5列的表格中，列的表格中，使其每行每列及对角线的和均为使其每行每列及对角线的和均为565。命令如下：命令如下：B=100+magic(5)B=100+magic(5)B=117 124 101 108 115 123 105 107 114 116 104 106 113 120 122 110 112 119 121 103 111 118 125 102 109 sum(B(1,:)ans=565 sum(B(2,:)ans=565 sum(B(:,4)ans=565 B(1,1)+B(2,2)+B(3,3)+B(4,4)+B(5,5)ans=5652.范得蒙矩阵范得蒙矩阵 函数函数vander(V)生成以向量生成以向量V为基础向量的范得为基础向量的范得蒙矩阵。蒙矩阵。VANDER Vandermonde matrix.A=VANDER(V)returns the Vandermonde matrix whose columns are powers of the vector V,that is A(i,j)=v(i)(n-j).p=1 2 3 4 5p=1 2 3 4 5 A=VANDER(p)A=1 1 1 1 1 16 8 4 2 1 81 27 9 3 1 256 64 16 4 1 625 125 25 5 1 3.希尔伯特矩阵 Hilbert矩阵的每个元素的值，由行数i和列数j决定，等于1/(i+j-1)，生成希尔伯特矩阵的函数是hilb(n)。MATLAB中，有一个专门求希尔伯特矩阵的逆的函数invhilb(n)，其功能是求n阶的希尔伯特矩阵的逆矩阵。hilb(4)ans=1.0000 0.5000 0.3333 0.2500 0.5000 0.3333 0.2500 0.2000 0.3333 0.2500 0.2000 0.1667 0.2500 0.2000 0.1667 0.1429 I=invhilb(4)I=16 -120 240 -140 -120 1200 -2700 1680 240 -2700 6480 -4200 -140 1680 -4200 28004.托普利兹矩阵托普利兹矩阵 生成托普利兹矩阵的函数是生成托普利兹矩阵的函数是toeplitz(x,y)，它生成一个以，它生成一个以x为第为第1列，列，y为第为第1行的托普利兹矩阵。这里行的托普利兹矩阵。这里x,y均为向量，均为向量，二者不必等长。当二者不必等长。当x和和y的第一个元素不同时，系统将给的第一个元素不同时，系统将给出提示信息并以出提示信息并以x中的元素为准中的元素为准 c=1 2 3 4 5;r=1.5 2.5 3.5 4.5 5.5;toeplitz(c,r)Warning:First element of input column does not match first element of input row.Column wins diagonal conflict.(Type warning off MATLAB:toeplitz:DiagonalConflict to suppress this warning.)In E:matlabanzhuangtoolboxmatlabelmattoeplitz.m at line 18ans=1.0000 2.5000 3.5000 4.5000 5.5000 2.0000 1.0000 2.5000 3.5000 4.5000 3.0000 2.0000 1.0000 2.5000 3.5000 4.0000 3.0000 2.0000 1.0000 2.5000 5.0000 4.0000 3.0000 2.0000 1.00005.友矩阵 Compan矩阵生成多项式系数向量P的伴随矩阵，其中（A(1,:)=-P(2:n)/P(1)）,友矩阵的函数是：compan(P)。P是一个多项式的系数向量，高次幂系数排在前，低次幂排在后。例5.3求x3-7x+6的根。u=1 0-7 6u=1 0 -7 6 A=compan(u)A=0 7 -6 1 0 0 0 1 0 eig(compan(u)ans=-3.0000 2.0000 1.0000该多项式的根是该多项式的根是-3,2,1友矩阵的特征根正好是多项友矩阵的特征根正好是多项式的根式的根6.帕斯卡矩阵帕斯卡矩阵 Pascal矩阵是一个实对称的正定矩阵，它由矩阵是一个实对称的正定矩阵，它由Pascal三角形组三角形组成，成，Pascal三角形是由三角形是由0到到2n-1阶的二项式系数组成，把二项阶的二项式系数组成，把二项式系数依次填写在矩阵的左侧对角线上，提取左侧的式系数依次填写在矩阵的左侧对角线上，提取左侧的n行行n 列即列即为为Pascal矩阵。矩阵。函数函数pascal(n)生成一个生成一个n阶的帕斯卡矩阵。阶的帕斯卡矩阵。阶数二项式系数0111 121 2 131 3 3 141 4 6 4 151 5 10 10 5 11 1 1 1 1 2 3 4 1 3 6 101 4 16 20 1 1 15 6151 5 151 61Pascal(4)例5.4求(x+y)5的展开式。在MATLAB命令窗口，输入命令：pascal(6)ans=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 252其其次对角线次对角线上的元素上的元素1，5，10，10，5，1即为展开即为展开式的系数。式的系数。5.2 矩阵分析5.2.1 矩阵结构变换1.矩阵的转置转置运算符是单撇号()。2.矩阵的旋转 矩阵的旋转利用函数rot90(A,k)，功能是将矩阵A旋转90的k倍，当k为1时可省略。3.矩阵的左右翻转对矩阵A实施左右翻转的函数是fliplr(A)。4.矩阵的上下翻转对矩阵A实施上下翻转的函数是flipud(A)。目录目录A=magic(3)A=8 1 6 3 5 7 4 9 2 fliplr(A)%将矩阵的列左右翻转将矩阵的列左右翻转ans=6 1 8 7 5 3 2 9 4 flipud(A)%将矩阵将矩阵A的行上、的行上、下翻转下翻转ans=4 9 2 3 5 7 8 1 6 rot90(A)%将矩阵将矩阵A逆时针旋转逆时针旋转90度度ans=6 7 2 1 5 9 8 3 4 5.2.2 矩阵的逆与伪逆1.矩阵的逆 求一个矩阵的逆非常容易。求方阵A的逆可调用函数inv(A)。例5.4 用求逆矩阵的方法解线性方程组。命令如下：一般情况下，用左除比求矩阵的逆的方法更有效，即x=Ab。目录目录 A=1,2,3;1,4,9;1,8,27;b=5,-2,6;x=inv(A)*bx=23.0000 -14.5000 3.6667 x=Abx=23.0000 -14.5000 3.6667 2.矩阵的伪逆MATLAB中，求一个矩阵伪逆的函数是pinv(A)。例5.5 求A的伪逆，并将结果送B。命令如下：A=3,1,1,1;1,3,1,1;1,1,3,1;B=pinv(A)例5.6 求矩阵A的伪逆。在MATLAB命令窗口，输入命令：A=0,0,0;0,1,0;0,0,1;pinv(A)目录目录B=0.3929 -0.1071 -0.1071 -0.1071 0.3929 -0.1071 -0.1071 -0.1071 0.3929 0.0357 0.0357 0.0357ans=0 0 0 0 1 0 0 0 1 5.2.3 方阵的行列式方阵的行列式求方阵求方阵A所对应的行列式的值的函数是所对应的行列式的值的函数是det(A)。例例5.7用克莱姆用克莱姆(Cramer)方法求解线性方程组。方法求解线性方程组。程序如下：程序如下：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;%定义系数矩阵定义系数矩阵b=4;6;12;6;%定义常数项向量定义常数项向量D1=b,D(:,2:4);%用方程组的右端向量置换用方程组的右端向量置换D的第的第1列列D2=D(:,1:1),b,D(:,3:4);%用方程组的右端向量置换用方程组的右端向量置换D的第的第2列列D3=D(:,1:2),b,D(:,4:4);%用方程组的右端向量置换用方程组的右端向量置换D的第的第3列列D4=D(:,1:3),b;%用方程组的右端向量置换用方程组的右端向量置换D的第的第4列列DD=det(D);x1=det(D1)/DD;x2=det(D2)/DD;x3=det(D3)/DD;x4=det(D4)/DD;x1,x2,x3,x4目录目录ans=1 1 -1 -1 5.2.4 矩阵的秩MATLAB中，求矩阵秩的函数是rank(A)。例如，求例5.7中方程组系数矩阵D的秩，命令是：说明D是一个满秩矩阵。目录目录 D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;r=rank(D)r=4 5.2.5 向量和矩阵的范数(不讲)1.计算向量3种常用范数的函数(1)norm(V)或norm(V,2)计算向量V的2范数=sum(abs(V).2)(1/2)(2)norm(V,1)=sum(abs(V)计算向量V的1范数(3)norm(V,inf)计算向量V的范数=max(abs(V)目录目录例5.8 已知V，求V的3种范数。v1=norm(V,1)%求V的1范数sum(abs(V)v1=5/2 v2=norm(V)%求V的2范数sum(abs(V).2)(1/2)v2=3/2 v3=norm(V,inf)%求V的范数max(abs(V)v3=1 2.矩阵的范数及其计算函数MATLAB中提供了求3种矩阵范数的函数，其函数调用格式与求向量的范数的函数完全相同NORM(X)is the largest singular value of X,max(svd(X).NORM(X,2)is the same as NORM(X).NORM(X,1)is the 1-norm of X,the largest column sum,=max(sum(abs(X).NORM(X,inf)is the infinity norm of X,the largest row sum,=max(sum(abs(X).A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19;a1=norm(A,1)%求A的1范数a2=norm(A)%求A的2范数ainf=norm(A,inf)%求A的范数例5.9 求矩阵A的三种范数。命令如下：a1=75 a2=2790/47 ainf=75 5.2.6 矩阵的条件数和迹1.的条件数MATLAB中，计算矩阵A的3种条件数的函数是：(1)cond(A,1)计算A的1范数下的条件数(2)cond(A)或cond(A,2)计算A的2范数数下的条件数(3)cond(A,inf)计算A的 范数下的条件数目录目录例5.10 求矩阵X的三种条件数。命令如下：A=2,2,3;4,5,-6;7,8,9;C1=cond(A,1)C2=cond(A)C3=cond(A,inf)C1=1044/7 C2=7126/81 C3=144 2.矩阵的迹MATLAB中，求矩阵的迹的函数是trace(A)。例如，X=2 2 3;4 5-6;7 8 9;trace(X)ans=16目录目录5.2.7 矩阵的特征值与特征向量MATLAB中，计算矩阵A的特征值和特征向量的函数是eig(A)，常用的调用格式有3种：(1)E=eig(A)求矩阵A的全部特征值，构成向量E。(2)V,D=eig(A)求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。(3)V,D=eig(A,nobalance)与第2种格式类似，但第2种格式中先对A作相似变换后求矩阵A的特征值和特征向量，而格式3直接求矩阵A的特征值和特征向量。例5.11 用3种不同的格式求A的特征值和特征向量。命令如下：A=1,2,2;1,-1,1;4,-12,1;E=eig(A)V,D=eig(A)V,D=eig(A,nobalance)目录目录例5.12用求特征值的方法解方程3x5-7x4+5x2+2x-18。命令如下：p=3,-7,0,5,2,-18;A=compan(p);%A的友矩阵x1=eig(A)%求A的特征值x2=roots(p)%直接求多项式p的零点两种方法求得的方程的根是完全一致的，实际上，roots函数正是应用求友矩阵的特征值的方法来求方程的根。x1=2.1837 1.0000+1.0000i 1.0000-1.0000i -0.9252+0.7197i -0.9252-0.7197ix2=2.1837 1.0000+1.0000i 1.0000-1.0000i -0.9252+0.7197i -0.9252-0.7197i5.2.8 MATLAB在三维向量中的应用1.向量共线或共面的判断例5.13 设X=(1，1，1)，Y=(-1，2，1)，Z=(2，2，2)，判断这三个向量的共线共面问题。命令如下：X=1,1,1;Y=-1,2,1;Z=2,2,2;XY=X;Y;YZ=Y;Z;ZX=Z;X;XYZ=X;Y;Z;rank(XY)rank(YZ)rank(ZX)rank(XYZ)目录目录ans=2ans=2ans=1ans=2 2.向量方向余弦的计算例5.14设向量V=(5，-3，2)，求V的方向余弦。建立一个函数文件direct.m：function f=f(v)r=norm(v);if r=0 f=0else f=v(1)/r,v(2)/r,v(3)/r;endreturn在MATLAB命令窗口，输入命令：v=5,-3,2;f=direct(v)目录目录f=2220/2737 -1332/2737 888/2737 3.向量的夹角例5.15 设U=(1，0，0)，V=(0，1，0)，求U,V间的夹角.命令如下：U=1,0,0;V=0,1,0;r1=norm(U);r2=norm(V);UV=U*V;cosd=UV/r1/r2;D=acos(cosd)4.两点间的距离例5.16 设 U=(1，0，0)，V=(0，1，0)，求U、V两点间的距离。命令如下：U=1,0,0;V=0,1,0;UV=U-V;D=norm(UV)目录目录 5.向量的向量积例5.17设U=(2，-3，1)，V=(3，0，4)，求UV。命令如下：U=2,-3,1;V=3,0,4;W=eye(3);A1=W(1,:);U;V;A2=W(2,:);U;V;A3=W(3,:);U;V;UV=det(A1),det(A2),det(A3)UV=-12 -5 96.向量的混合积例5.18 设U=(0，0，2)，V=(3，0，5)，W=(1，1，0)，求以这三个向量构成的六面体的体积。命令如下：U=0,0,2;V=3,0,5;W=1,1,0;A=U;V;W;det(A)ans=6目录目录7.点到平面的距离例5.19求原点到平面X+Y+Z=1的距离。命令如下：u=0,0,0;v=1,1,1;%A=B=C=1,u1=u2=u3=0,D=-1r=abs(u*v-1)/norm(v,2)r=0.5774目录目录5.3 矩阵分解与线性方程组求解矩阵分解与线性方程组求解5.3.1矩阵分解1.实对称矩阵的QDQ分解例5.20设对称矩阵A，对A进行QDQ分解。命令如下：A=2,1,4,6;1,2,1,5;4,1,3,4;6,5,4,2;Q,D=eig(A)Q*D*Qans=2.0000 1.0000 4.0000 6.0000 1.0000 2.0000 1.0000 5.0000 4.0000 1.0000 3.0000 4.0000 6.0000 5.0000 4.0000 2.0000结果与A相等，说明确实将A分解为了QDQ的乘积。目录目录 例5.21求下列二次型的标准形式及变换矩阵。命令如下：A=1,2,1;2,1,1;1,1,3;Q,D=eig(A)进一步作线性变换即得关于u,v,w的标准二次型：2.矩阵的LU分解MATLAB中，完成LU分解的函数是：(1)L,U=lu(A)将方阵A分解为交换下三角矩阵L和上三角矩阵U，使 A=LU。(2)L,U,P=lu(A)将方阵A分解为下三角矩阵L和上三角矩阵U，使 PA=LU。LU分解常用于求行列式以及解线性方程组。目录目录3.矩阵的QR分解（正交分解）对矩阵A进行QR分解的函数是Q,R=qr(A)，根据方阵A，求一个正交矩阵Q和一个上三角矩阵R，使A=Q*R。例如，对矩阵A进行QR分解的命令是：A=2,1,-2;1,2,1;2,5,3;Q,R=qr(A)目录目录 5.3.2 线性方程组求解线性方程组求解1.线性方程组解的一般讨论线性方程组解的一般讨论解线性方程组的一般解线性方程组的一般函数文件函数文件如下：如下：function x,y=line_solution(A,b)m,n=size(A);y=;if norm(b)0%非齐次方程组 if rank(A)=rank(A,b)%方程组相容 if rank(A)=m%有唯一解 x=Ab;else%方程组有无穷多个解,基础解系 disp(原方程组有有无穷个解，其齐次方程组的基础解系为y，特解为x);y=null(A,r);x=Ab;end else%方程组不相容，给出最小二乘法解 disp(方程组的最小二乘法解是：);x=Ab;目录目录end else%齐次方程组 if rank(A)=n%列满秩 x=zero(m,1)%0解 else%非0解 disp(方程组有无穷个解，基础解系为x);x=null(A,r);end endreturn2.应用举例例5.23求线性方程组的解。在MATLAB命令窗口，输入命令：A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6;x,y=line_solution(A,b)%调用自定义函数目录目录x=1 1 -1 -1 y=例例5.24求下列线性方程组的解。在MATLAB命令窗口，输入命令：A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,2;x,y=line_solution(A,b)原方程组有有无穷个解，其齐次方程组的基础解系为y，特解为xWarning:Rank deficient,rank=2 tol=8.6112e-015.In E:matlabanzhuangworkline_solution.m at line 10 x=-2/11 10/11 0 0 y=1/11 -9/11 -5/11 1/11 1 0 0 1 5.4 数据处理与多项式计算数据处理与多项式计算5.4.1 数据统计与分析1.求矩阵最大和最小元素(1)求向量的最大最小元素y=max(X)返回向量X的最大元素存入y。y,I=max(X)返回向量X的最大元素存入y，最大元素的序号存入I。(2)求矩阵的最大和最小元素max(A)返回一个行向量，向量的第i个元素是A矩阵的第i列上的最大元素。Y,U=max(A)返回两个行向量，Y向量记录A的每列的最大元素，U向量记录每列最大元素的行号。max(A,dim)dim取1或2。dim取1时，该函数和max(A)完全相同。dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大元素。目录目录(3)两个向量或矩阵对应元素的比较U=max(A,B)A,B是两个同型的向量或矩阵。结果U是与A,B同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。U=max(A,n)n是一个标量。结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。min函数的用法和max完全相同。目录目录例5.25 求矩阵A的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元。命令如下：A=13,-56,78;25,63,-235;78,25,563;1,0,-1;max(A,2)%求每行最大元素min(A,2)%求每行最小元素max(A)%求每列最大元素min(A)%求每列最小元素max(max(A)%求整个矩阵的最大元素min(min(A)%求整个矩阵的最小元素目录目录ans=78 63 563 1 ans=-56 -235 25 -1 ans=78 63 563 ans=1 -56 -235 ans=563 ans=-235 2.求矩阵的平均值和中值 求矩阵和向量元素的平均值的函数是mean，求中值的函数是median。它们的调用方法和max函数完全相同。3.矩阵元素求和与求积矩阵和向量求和与求积的基本函数是sum和prod，其使用方法和max类似。目录目录例5.26求矩阵A的每行元素的乘积和全部元素的乘积。命令如下：A=1,2,3,4;5,6,7,8;9,10,11,12;S=prod(A,2)prod(S)%求A的全部元素的乘积目录目录4.矩阵元素累加和与累乘积MATLAB中，使用cumsum和cumprod函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的用法和sum及prod相同例5.27求向量X=(1！，2！，3！，10！)。命令如下：Format long;X=cumprod(1:10)X=1 2 6 24 120 720 5040 40320 362880 36288005.标准方差 MATLAB中，提供了计算数据序列的标准方差的函数std。对于向量X，std(X)返回一个标准方差。对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵A各列或各行的标准方差。std函数的一般调用格式为：std(A,FLAG,dim)其中dim取1或2。当dim=1时，求各列元素的标准方差；当dim=2时，则求各行元素的标准方差。FLAG取0或1。目录目录6.元素排序 MATLAB中对向量X是排序函数是sort(X)，函数返回一个对X中的元素按升序排列的新向量。sort函数也可以对矩阵A的各列(或行)重新排序，其调用格式为：Y,I=sort(A,dim)其中dim指明对A的列还是行进行排序，若dim=1，则按列排，若dim=2，则按行排。Y是排序后的矩阵，而I记录Y中的元素在A中位置。目录目录例5.28对矩阵做各种排序。命令如下：A=1,-8,5;4,12,6;13,7,-13;sort(A)%对A的每列按升序排序-sort(-A,2)%对A的每行按降序排序X,I=sort(A)%对A按列排序，并将每个元素所在行号送矩阵I目录目录ans=1 -8 -13 4 7 5 13 12 6ans=5 1 -8 12 6 4 13 7 -13X=1 -8 -13 4 7 5 13 12 6I=1 1 3 2 3 1 3 2 25.4.2 数值插值数值插值1.一维数值插值一维数值插值 interp1函数调用格式为：Y1=interp1(X,Y,X1,method)函数根据X、Y的值，计算函数在X1处的值。X、Y是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，Y1是一个与X1等长的插值结果。method是插值方法，允许的取值有linear(线性插值)、nearest(最近插值)、spline(三次样条插值)、cubic（三次多项式插值），缺省值是linear。目录目录例例5.29用不同的插值方法计算sin(x)在/2点的值。这是一个一维插值问题。在MATLAB命令窗口，输入命令：X=0:0.2:pi;Y=sin(X);%给出X、Yinterp1(X,Y,pi/2)%用缺省方法(即线性插值方法)计算sin(/2)interp1(X,Y,pi/2,nearest)%用最近方法计算sin(/2)interp1(X,Y,pi/2,linear)%用线性方法计算sin(/2)interp1(X,Y,pi/2,spline)%用三次样条方法计算sin(/2)interp1(X,Y,pi/2,cubic)%用三次多项式方法计算sin(/2)MATLAB中有一个专门的三次样条插值函数Y1=spline(X,Y,X1)，其功能及使用方法与函数Y1=interp1(X,Y,X1,spline)完全相同。目录目录 例5.30 已知检测参数f随时间t的采样结果，用数值插值法计算t=2，7，12，17，22，17，32，37，42，47，52，57时f的值。这是一个一维数值插值问题，命令如下：T=0:5:65;X=2:5:57;F=3.2019,2.2560,879.5,1835.9,2968.8,4136.2,5237.9,6152.7,6725.3,6848.3,6403.5,6824.7,7328.5,7857.6;F1=interp1(T,F,X)%用线性方法插值F1=interp1(T,F,X,nearest)%用最近方法插值F1=interp1(T,F,X,spline)%用三次样条方法插值 F1=interp1(T,F,X,cubic)%用三次多项式方法插值目录目录2.二维数值插值 MATLAB中，提供了解决二维插值问题的函数。其调用格式为：Z1=interp2(X,Y,Z,X1,Y1,method)其中X、Y是两个向量，分别描述两个参数的采样点，Z是与参数采样点对应的采样变量的样本值，X1、Y1是两个向量或标量，描述欲插值的点。method的取值与一维插值函数相同。目录目录例5.31设Z=x2+y2，对Z函数在(0，1)(0，2)区域内进行插值。命令如下：x=0:0.1:10;y=0:0.2:20;X,Y=meshgrid(x,y);Z=X.2+Y.2;interp2(x,y,Z,0.5,0.5)%对函数在(0.5,0.5)点进行插值interp2(x,y,Z,0.5 0.6,0.4)%对函数在(0.5,0.4)点和(0.6，0.4)点进行插值interp2(x,y,Z,0.5 0.6,0.4 0.5)%对函数在(0.5,0.4)点和(0.6，0.5)点进行插值interp2(x,y,Z,0.5 0.6,0.4 0.5)%对函数在(0.5,0.4),(0.6,0.4),(0.5,0.5)和(0.6,0.5)点进行插值目录目录ans=0.5100ans=0.4100 0.5200ans=0.4100 0.6200ans=0.4100 0.5200 0.5100 0.62003.三维数值插值 对三维函数插值的函数是interp3，其使用方法和interp2相同。其调用格式为：W1=interp3(X,Y,Z,W,X1,Y1,Z1,method)函数返回三维插值结果。其中X、Y、Z是三个向量，分别描述三个参数的采样点，W是与参数采样点对应的采样变量的样本值，X1、Y1、Z1是三个向量或标量，描述欲插值的点。method是插值方法，可选，其缺省值是 line。method的取值与一、二维插值函数相同。目录目录5.4.3 曲线拟合 MATLAB中，提供了解决使用最小二乘法进行曲线拟合的函数。调用格式为：P,S=polyfit(X,Y,m)函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X、Y是两个等长的向量，P是一个长度为m+1的向量。目录目录例5.32 用一个5次多项式在区间0，2内逼近函数sin(x)。命令如下：目录目录 X=linspace(0,2*pi,50);Y=sin(X);P,S=polyfit(X,Y,5)%得到5次多项式的系数和误差plot(X,Y,k*,X,polyval(P,X),k-)P=-65/11683 287/3284 -58/147 578/2153 1528/1737 143/13993 S=R:6x6 double df:44 normr:132/3913 5.4.4 多项式计算1.多项式的建立 已知一个多项式的全部根X求多项式系数的函数是poly(X)，该函数返回以X为全部根的一个多项式P，当X是一个长度为m的向量时，P是一个长度为m+1的向量。目录目录 R=-1-2-3-4+3*j-4-3*j;%输入多项式根的向量输入多项式根的向量 P=poly(R)%构造多项式的系数向量构造多项式的系数向量P=1 14 84 244 323 1502.多项式求根多项式求根求多项式求多项式p(x)的根的函数是的根的函数是roots(P)，这里，这里，P是是p(x)的系数向量，该函数返回方程的系数向量，该函数返回方程p(x)=0的全部根的全部根(含含重根，复根重根，复根)。例：已知例：已知5阶多项式的系数向量为阶多项式的系数向量为P1828586730，求多项式的根，求多项式的根 P=1 8 28 58 67 30;R=roots(P)R=-1.00000000000000+2.00000000000000i-1.00000000000000-2.00000000000000i-2.99999999999998 -2.00000000000001 -1.00000000000000 3.多项式求值多项式求值求多项式求多项式p(x)在某点或某些点的函数值的函数是在某点或某些点的函数值的函数是polyval(P,x)。若。若x为一数值，则求多项式在该点的为一数值，则求多项式在该点的值；若值；若x为向量或矩阵，则对向量或矩阵中的每个元为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。素求其多项式的值。例例5.33 已知一个多项式已知一个多项式3x5+4x3-5x2-7.2x+5，计算：，计算：(1)计算计算f(x)=0 的全部根。的全部根。(2)由方程由方程f(x)=0的根构造一个多项式的根构造一个多项式g(x)，并与，并与f(x)进行对比。进行对比。(3)计算计算f(5)、f(7.8)、f(9.6)、f(12.3)的值。的值。命令如下：命令如下：P=3,0,4,-5,-7.2,5;X=roots(P)%求方程求方程f(x)=0的根的根G=poly(X)%求多项式求多项式g(x)X0=5,7.8,9.6,12.3;f=polyval(P,X0)%求多项式求多项式f(x)在给定点的值在给定点的值 多项式求值还有一个函数是多项式求值还有一个函数是polyvalm，其调用格式与，其调用格式与polyval相相同，但含义不同。同，但含义不同。polyvalm函数要求函数要求x为方阵，它以方阵为自变量求为方阵，它以方阵为自变量求多项式的值。多项式的值。目录目录X=-0.30455728045903+1.62172723875407i-0.30455728045903-1.62172723875407i-1.00664150535070 1.01901848386455 0.59673758240420 G=1.00000000000000 0.00000000000000 1.33333333333333 -1.66666666666667 -2.39999999999999 1.66666666666666f=1.0e+005*0.09719000000000 0.88158079040000 2.47625833280000 8.51195163290000 4.多项式的四则运算(1)多项式的加减法(2)多项式的乘法函数conv(P1,P2)用于求多项式P1和P2的乘积。(3)多项式的除法函数Q,r=deconv(P1,P2)用于对多项式P1和P2作除法运算。其中Q返回多项式P1除以P2的商式，r返回P1除以P2的余式。这里，Q和r仍是多项式系数向量。deconv是conv的逆函数，即有P1=conv(P2,Q)+r。目录目录例5.34设有两个多项式，计算：(1)求f(x)+g(x)、f(x)-g(x)。(2)求f(x)g(x)、f(x)/g(x)。目录目录在在MATLAB命令窗口，输入命令命令窗口，输入命令：f=3,-5,2,-7,5,6;g=3,5,-3;g1=0,0,0,g;f+g1%求f(x)+g(x)f-g1%求f(x)-g(x)conv(f,g)%求f(x)*g(x)Q,r=deconv(f,g)%求f(x)/g(x)，商式送Q，余式送r。5.多项式的导函数对多项式求导数的函数是：p=polyder(P)求多项式P的导函数p=polyder(P,Q)求P*Q的导函数p,q=polyder(P,Q)求P/Q的导函数，导函数的分子存入p，分母存入q。目录目录例5.35求有理分式 的导数。命令如下：目录目录P=3,5,0,-8,1,-5;Q=10,5,0,0,6,0,0,7,-1,0,-100;p,q=polyder(P,Q)5.4.5 函数的最大值与最小值MATLAB中用于求最小值的函数是：fmin(f,a,b)求单变量函数f(x)在区间(a,b)上的最小值点。fmins(F,X0)求多变量函数F(x)在估计值X0附近的最小值点。MATLAB没有专门提供求函数最大值点的函数，但只要注意到-f(x)在区间(a,b)上的最小值点就是f(x)在(a,b)的最大值点，所以fmin(-f,a,b)返回函数f(x)在区间(a,b)上的最大值。目录目录例5.36 求函数f(x)在区间(-10,-1)和(1,10)上的最小值点。首先建立函数文件fx.m：function f=f(x)f=x-1/x+5;return再在MATLAB命令窗口，输入命令：fmin(fx,-10,-1)%求函数在区间(-10，-1)内的最小值点fmin(f,1,10)%求函数在区间(1，10)内的最小值点。注意函数名f不用加目录目录例5.37 设有函数f(x,y,z)，求函数f在(0.5,0.5,0.5)附近的最小值。建立函数文件fxyz.m：function f=f(u)x=u(1);y=u(2);z=u(3);f=x+y.2./x/4+z.2./y+2./z;return在MALAB命令窗口，输入命令：U=fmins(fxyz,0.5,0.5,0.5)%求函数的最小值点fxyz(U)%求函数的最小值目录目录5.5 傅立叶分析(不讲)MATLAB中，提供了对向量(或直接对矩阵的行或列)进行离散傅立叶变换的函数，其调用格式是：Y=fft(X,n,dim)(1)当X是一个向量时，返回对X的离散傅立叶变换。(2)当X是一个矩阵时，返回一个矩阵并送Y，其列(行)是对X的列(行)的离散傅立叶变换。目录目录例5.38 求X=(1，0，-3，5，2)的离散傅立叶逆变换。在MATLAB命令窗口，输入命令：X=1,0,-3,5,2;Y=fft(X)%对X进行变换3.离散傅立叶变换的逆变换MATLAB中，对向量(或直接对矩阵的行或列)进行离散傅立叶逆变换的函数的调用方法是：Y=ifft(X,n,dim)函数对X进行离散傅立叶逆变换。其中X、n、dim的意义及用法和离散傅立叶变换函数fft完全相同。目录目录例5.39 对矩阵A的列向量、行向量分别进行离散傅立叶变换、并对变换结果进行逆变换。命令如下：A=3,2,1,1;-5,1,0,1;3,2,1,5;fftA=fft(A)%求A的列向量的傅立叶变换fftA2=fft(A,4,2)%求A的行向量的傅立叶变换ifft(fftA)%对矩阵fftA的列向量进行傅立叶逆变换，结果应等于Aifft(fftA2,4,2)%对矩阵fftA2的行向量进行傅立叶逆变换，其结果应等于A目录目录5.6 数值微积分5.6.1 数值微分MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数。DX=diff(X)计算向量X的向前差分，DX(i)=X(i+1)-X(i)，0in。DX=diff(X,n)计算X的n阶向前差分，diff(X,2)=diff(diff(X)。DX=diff(A,n,dim)计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分，dim=2，按行计算差分。目录目录例5.40 求向量sin(X)的13阶差分。设X由0，2间均匀分布的10个点组成。命令如下：X=linspace(0,2*pi,10);Y=sin(X);DY=diff(Y);%计算Y的一阶差分D2Y=diff(Y,2);%计算Y的二阶差分，也可用命令diff(DY)计算D3Y=diff(Y,3);%计算Y的三阶差分，也可用diff(D2Y)或diff(DY,2)目录目录例5.41 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f(x)的图象。程序如下：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5);%用5次多项式p拟合f(x)dp=polyder(p);%对拟合多项式p求导数dpdpx=polyval(dp,x);%求dp在假设点的函数值dx=diff(f(x,3.01)/0.01;%直接对f(x)求数值导数gx=g(x);%求函数f的导函数g在假设点的导数plot(x,dpx,x,dx,g.,x,gx,r-);%作图目录目录5.6.2数值积分(1)被积函数是一个解析式 函数quad(f,a,b,tol,trace)用于求被积函数f(x)在a,b上的定积分，tol是计算精度，缺省值是0.001。trace非0时，画出积分图形。注意，调用quad函数时，先要建立一个描述被积函数f(x)的函数文件或语句函数。当被积函数f含有一个以上的变量时，quad函数的调用格式为：quad(f,a,b,tol,trace,g1,g2)其中f,a,b,tol,trace等参数的含义同前。数值积分函数还有一种形式quad8，其用法与quad完全相同。目录目录例5.42 用两种不同的方法求积分。先建立一个函数文件ex.m：function ex=ex(x)ex=exp(-x.2);%注意应用点运算return然后，在MATLAB命令窗口，输入命令：quad(ex,0,1,1e-6)%注意函数名应加字符引号quad8(ex,0,1,1e-6)%用另一函数求积分目录目录(2)被积函数由一个表格定义 MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X、Y定义函数关系Y=f(X)。例5.43用trapz函数计算积分。在MATLAB命令窗口，输入命令：X=0:0.01:1;Y=exp(-X.2);trapz(X,Y)目录目录(3)二重积分例5.44计算二重积分。建立一个函数文件fixy.m：function f=f(x,y)f=exp(-x.2-y.2);return建立一个命令文件ftxy1.m：for i=1:20 int2(i)=quad(fixy,0,1,x(i);%在二维函数fixy中以x=x(i)代入并对y积分。end在MATLAB命令窗口，输入命令：x=linspace(0,1,20);ftxy1trapz(x,int2)目录目录实际上，MATLAB提供了计算二重积分的函数：dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在a,bc,d区域上的二重积分。参数tol，trace的用法与函数quad完全相同。如果直接使用这里介绍的二重积分函数dblquad来求解本例就非常简单，命令如下：g=inline(exp(-x.2-y.2);dblquad(g,0,1,0,1)%直接调用二重积分函数求解目录目录5.7 常微分方程的数值求解基于龙格库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为：X,Y=ode23(f,x0,xn,y0)X,Y=ode45(f,x0,xn,y0)其中X、Y是两个向量，X对应自变量x在求解区间x1，xn的一组采样点，其采样密度是自适应的，无需指定；Y是与X对应的一组解，f是一个函数，x0，xn代表自变量的求解区间，y0=y(x0)，由方程的初值给定。函数在求解区间x0，xn内，自动设立采样点向量X，并求出解函数y在采样点X处的样本值。目录目录例5.45 求微分方程初值问题在1，3区间内的数值解，并将结果与解析解进行比较。先建立一个该函数的m文件fxy1.m：function f=f(x,y)f=-2.*y./x+4*x%注意使用点运算符return再输入命令：X,Y=ode45(fxy1,1,3,2);X%显示自变量的一组采样点Y%显示求解函数与采样点对应的一组数值解(X.2+1./X.2)%显示求解函数与采样点对应的一组解析解目录目录例5.46 求解初值问题在区间0，2中的解。建立一个函数文件 fxy2.m：function f=f(x,y)f(2)=-x.*y(2)+x.2-5;f(1)=y(2);f=f;return在MATLAB命令窗口，输入命令：X,Y=ode45(fxy2,0,2,5,6);X,Y目录目录5.8 非线性方程的数值求解1单变量非线性方程求解MATLAB中，提供了求解单变量方程的函数fzero(f,x0,tol)，该函数采用迭代法计算函数f(x)的一个零点，迭代初值为x0，当两次迭代结果小于tol时停止迭代过程。tol的缺省值是eps。注意，在调用函数fzero 之前，要使用m文件建立自己要计算的函数f(x)，只有定义了函数f(x)的m文件后，才能在fzero函数的参数中使用自定义函数名。目录目录例5.47 求f(x)=x-+5 在x0=-5和x0=1作为迭代初值时的零点。先编制一个函数文件fz.m：function f=f(x)f=x-1/x+5;然后，在MATLAB命令窗口，输入命令：fzero(fz,-5)%以-5作为迭代初值Zero found in the interval:-4.8,-5.2.fzero(fz,1)目录目录2非线性方程组求解函数fsolve调用格式为：X=fsolve(F,X0)例5.48 求方程组在(1，1，1)附近的解并对结果进行验证。首先建立方程的函数文件fxyz1.m：function F=F(X)x=X(1);y=X(2);z=X(3);F(1)=sin(x)+y+z2*exp(x);F(2)=x+y*z;F(3)=x*y*z;在MATLAB命令窗口，输入命令：X=fsolve(fxyz1,1,1,1)%求解X的三个分量x、y、zY=fxyz1(X)%检验所求结果X是否满足原方程组norm(Y)%求Y向量的模目录目录例5.49 求圆和直线的两个交点。建立方程组函数文件fxyz2.m：function F=F(X)x=X(1);y=X(2);z=X(3);F(1)=x2+y2+z2-9;F(2)=3*x+5*y+6*z;F(3)=x-3*y-6*z-1