《MATLAB在高等数学中的应用》由会员分享，可在线阅读，更多相关《MATLAB在高等数学中的应用（83页珍藏版）》请在金锄头文库上搜索。
1、第三章 MATLAB在高等数学中的应用电子信息学院8/2/2024电子信息学院3.1矩阵分析矩阵分析3.1.1 对角阵与三角阵对角阵与三角阵1对角阵对角阵只有对角线上有非只有对角线上有非0元素的矩阵称为对角元素的矩阵称为对角矩阵，对角线上的元素相等的对角矩阵矩阵，对角线上的元素相等的对角矩阵称为数量矩阵，对角线上的元素都为称为数量矩阵，对角线上的元素都为1的的对角矩阵称为单位矩阵。对角矩阵称为单位矩阵。8/2/2024电子信息学院(1) 提取矩阵的对角线元素提取矩阵的对角线元素设设A为为mn矩阵，矩阵，diag(A)函数用于提取矩阵函数用于提取矩阵A主对角线元主对角线元素，产生一个具有素，产生
2、一个具有min(m,n)个元素的列向量。个元素的列向量。diag(A)函数还有一种形式函数还有一种形式diag(A,k)，其功能是提取第，其功能是提取第k条条对角线的元素。对角线的元素。(2) 构造对角矩阵构造对角矩阵设设V为具有为具有m个元素的向量，个元素的向量，diag(V)将产生一个将产生一个mm对角对角矩阵，其主对角线元素即为向量矩阵，其主对角线元素即为向量V的元素。的元素。diag(V)函数也有另一种形式函数也有另一种形式diag(V,k)，其功能是产生一个，其功能是产生一个nn(n=m+)对角阵，其第对角阵，其第k条对角线的元素即为向量条对角线的元素即为向量V的元的元素。素。8/2
3、/2024电子信息学院例例 先建立先建立55矩阵矩阵A，然后将，然后将A的第一行元的第一行元素乘以素乘以1，第二行乘以，第二行乘以2，第五行乘，第五行乘以以5。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:5);D*A %用用D左乘左乘A，对，对A的每行的每行乘以一个指定常数乘以一个指定常数8/2/2024电子信息学院2三角阵三角阵三角阵又进一步分为上三角阵和下三角三角阵又进一步分为上三角阵和下三角阵，所谓上三角阵，即矩阵的对角线以阵，所谓上三角阵，即矩阵的对角线以下的元素全为下的元素全
4、为0的一种矩阵，而下三角阵的一种矩阵，而下三角阵则是对角线以上的元素全为则是对角线以上的元素全为0的一种矩阵。的一种矩阵。8/2/2024电子信息学院(1)上三角矩阵上三角矩阵求矩阵求矩阵A的上三角阵的的上三角阵的MATLAB函数是函数是triu(A)。triu(A)函数也有另一种形式函数也有另一种形式triu(A,k)，其功，其功能是求矩阵能是求矩阵A的第的第k条对角线以上的元素。例条对角线以上的元素。例如，提取矩阵如，提取矩阵A的第的第2条对角线以上的元素，条对角线以上的元素，形成新的矩阵形成新的矩阵B。(2)下三角矩阵下三角矩阵在在MATLAB中，提取矩阵中，提取矩阵A的下三角矩阵的的下
5、三角矩阵的函数是函数是tril(A)和和tril(A,k)，其用法与提取上三，其用法与提取上三角矩阵的函数角矩阵的函数triu(A)和和triu(A,k)完全相同。完全相同。8/2/2024电子信息学院3.1.2 矩阵的转置与旋转矩阵的转置与旋转1矩阵的转置矩阵的转置转置运算符是单撇号转置运算符是单撇号()。2矩阵的旋转矩阵的旋转利用函数利用函数rot90(A,k)将矩阵将矩阵A旋转旋转90的的k倍，当倍，当k为为1时可省略。时可省略。8/2/2024电子信息学院3矩阵的左右翻转矩阵的左右翻转对矩阵实施左右翻转是将原矩阵的第一对矩阵实施左右翻转是将原矩阵的第一列和最后一列调换，第二列和倒数第二
6、列和最后一列调换，第二列和倒数第二列调换，列调换，依次类推。，依次类推。MATLAB对矩对矩阵阵A实施左右翻转的函数是实施左右翻转的函数是fliplr(A)。4矩阵的上下翻转矩阵的上下翻转MATLAB对矩阵对矩阵A实施上下翻转的函数实施上下翻转的函数是是flipud(A)。8/2/2024电子信息学院3.1.3 矩阵的逆与伪逆矩阵的逆与伪逆1矩阵的逆矩阵的逆对于一个方阵对于一个方阵A，如果存在一个与其同阶的方阵，如果存在一个与其同阶的方阵B，使，使得：得：AB=BA=I (I为单位矩阵为单位矩阵)则称则称B为为A的逆矩阵，当然，的逆矩阵，当然，A也是也是B的逆矩阵。的逆矩阵。求一个矩阵的逆是一
7、件非常烦琐的工作，容易出错，求一个矩阵的逆是一件非常烦琐的工作，容易出错，但在但在MATLAB中，求一个矩阵的逆非常容易。求方阵中，求一个矩阵的逆非常容易。求方阵A的逆矩阵可调用函数的逆矩阵可调用函数inv(A)。例例 用求逆矩阵的方法解线性方程组。用求逆矩阵的方法解线性方程组。Ax=b其解为：其解为：x=A-1b8/2/2024电子信息学院2矩阵的伪逆矩阵的伪逆如果矩阵如果矩阵A不是一个方阵时，矩阵不是一个方阵时，矩阵A没有没有逆矩阵，但可以找到一个与逆矩阵，但可以找到一个与A的转置矩阵的转置矩阵A同型的矩阵同型的矩阵B，使得：，使得：ABA=ABAB=B此时称矩阵此时称矩阵B为矩阵为矩阵A
8、的伪逆，也称为广的伪逆，也称为广义逆矩阵。在义逆矩阵。在MATLAB中，求一个矩阵中，求一个矩阵伪逆的函数是伪逆的函数是pinv(A)。8/2/2024电子信息学院3.1.4 方阵的行列式方阵的行列式把一个方阵看作一个行列式，并对其按把一个方阵看作一个行列式，并对其按行列式的规则求值，这个值就称为矩阵行列式的规则求值，这个值就称为矩阵所对应的行列式的值。在所对应的行列式的值。在MATLAB中，中，求方阵求方阵A所对应的行列式的值的函数是所对应的行列式的值的函数是det(A)。8/2/2024电子信息学院3.1.5 线性方程组求解线性方程组求解3.1.5.1 直接解法直接解法1利用左除运算符的直
9、接解法利用左除运算符的直接解法对于线性方程组对于线性方程组Ax=b，可以利用左除运算符，可以利用左除运算符“”求解：求解： x=Ab例例 用直接解法求解下列线性方程组。用直接解法求解下列线性方程组。命令如下：命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;x=Ab8/2/2024电子信息学院3.1.5.2利用矩阵的分解求解线性方程组利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有常见的矩阵分解有L
10、U分解、分解、QR分解、分解、Cholesky分解，以及分解，以及Schur分解、分解、Hessenberg分解、奇异分解等。分解、奇异分解等。8/2/2024电子信息学院(1) LU分解分解矩阵的矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵明，只要方阵A是非奇异的，是非奇异的，LU分解总是可以进行的。分解总是可以进行的。MATLAB提供的提供的lu函数用于对矩阵进行函数用于对矩阵进行LU分解，其调用分解，其调用格式为：格式为：L,U=lu(X)
11、：产生一个上三角阵：产生一个上三角阵U和一个变换形式的下和一个变换形式的下三角阵三角阵L(行交换行交换)，使之满足，使之满足X=LU。注意，这里的矩。注意，这里的矩阵阵X必须是方阵。必须是方阵。L,U,P=lu(X)：产生一个上三角阵：产生一个上三角阵U和一个下三角阵和一个下三角阵L以以及一个置换矩阵及一个置换矩阵P，使之满足，使之满足PX=LU。当然矩阵。当然矩阵X同样同样必须是方阵。必须是方阵。实现实现LU分解后，线性方程组分解后，线性方程组Ax=b的解的解x=U(Lb)或或x=U(LPb)，这样可以大大提高运算速度。，这样可以大大提高运算速度。8/2/2024电子信息学院例用例用LU分解
12、求解分解求解p79例例3-5线性方程组。线性方程组。A=6,3,4;-2,5,7;8,-4,-3;b=3,-4,-7;L,U=lu(A);x=U(Lb)或采用或采用LU分解的第分解的第2种格式，命令如下：种格式，命令如下：L,U ,P=lu(A);x=U(LP*b)8/2/2024电子信息学院 (2) QR分解分解对矩阵对矩阵X进行进行QR分解，就是把分解，就是把X分解为一个正交矩阵分解为一个正交矩阵Q和和一个上三角矩阵一个上三角矩阵R的乘积形式。的乘积形式。QR分解只能对方阵进分解只能对方阵进行。行。MATLAB的函数的函数qr可用于对矩阵进行可用于对矩阵进行QR分解，分解，其调用格式为：其
13、调用格式为：Q,R=qr(X)：产生一个一个正交矩阵：产生一个一个正交矩阵Q和一个上三角矩和一个上三角矩阵阵R，使之满足，使之满足X=QR。Q,R,E=qr(X)：产生一个一个正交矩阵：产生一个一个正交矩阵Q、一个上三角、一个上三角矩阵矩阵R以及一个置换矩阵以及一个置换矩阵E，使之满足，使之满足XE=QR。实现实现QR分解后，线性方程组分解后，线性方程组Ax=b的解的解x=R(Qb)或或x=E(R(Qb)。8/2/2024电子信息学院例例 用用QR分解求解线性方程组。分解求解线性方程组。命令如下：命令如下：A=6,3,4;-2,5,7;8,-4,-3;b=3,-4,-7;Q,R=qr(A);x
14、=R(Qb)或采用或采用QR分解的第分解的第2种格式，命令如下：种格式，命令如下：Q,R,E=qr(A);x=E*(R(Qb)8/2/2024电子信息学院3.1.5.3迭代解法迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括分析中，迭代解法主要包括 Jacobi迭代法、迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。迭代法、超松弛迭代法和两步迭代法。1Jacobi迭代法迭代法对于线性方程组对于线性方程组Ax=b，如果，如果A为非奇异方阵，则可将为非奇异方阵，则可将A分分解为解为A=D-L-U，其中
15、，其中D为对角阵，其元素为为对角阵，其元素为A的对角元的对角元素，素，L与与U为为A的下三角阵和上三角阵，于是的下三角阵和上三角阵，于是Ax=b化为：化为：x=D-1(L+U)x+D-1b与之对应的迭代公式为：与之对应的迭代公式为：x(k+1)=D-1(L+U)x(k)+D-1b这就是这就是Jacobi迭代公式。如果序列迭代公式。如果序列x(k+1)收敛于收敛于x，则，则x必是方程必是方程Ax=b的解。的解。8/2/2024电子信息学院Jacobi迭代法的迭代法的MATLAB函数文件函数文件Jacobi.m如下：如下：function y,n=jacobi(A,b,x0,eps)if narg
16、in=3 eps=1.0e-6;elseif nargin=eps x0=y; y=B*x0+f; n=n+1;end8/2/2024电子信息学院例例 用用Jacobi迭代法求解线性方程组。设迭代初值为迭代法求解线性方程组。设迭代初值为0，迭，迭代精度为代精度为10-6。在命令中调用函数文件在命令中调用函数文件Jacobi.m，命令如下：，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=jacobi(A,b,0,0,0,1.0e-6)(设设x1,x2,x3为为0,n为迭代的次数为迭代的次数)8/2/2024电子信息学院2Gauss-Serdel迭代法迭代
17、法将在将在Jacobi迭代过程中迭代过程中,原来的迭代公式原来的迭代公式Dx(k+1)=(L+U)x(k)+b改进为改进为 Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：，于是得到： x(k+1)=(D-L)-1Ux(k)+(D-L)-1b该式即为该式即为Gauss-Serdel迭代公式。和迭代公式。和Jacobi迭代相比，迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。迭代用新分量代替旧分量，精度会高些。8/2/2024电子信息学院Gauss-Serdel迭代法的迭代法的MATLAB函数文件函数文件gauseidel.m如下：如下：function y,n=
18、gauseidel(A,b,x0,eps)if nargin=3 eps=1.0e-6;elseif nargin=eps x0=y; y=G*x0+f; n=n+1;end8/2/2024电子信息学院例例 用用Gauss-Serdel迭代法求解下列线性方程组。设迭代迭代法求解下列线性方程组。设迭代初值为初值为0，迭代精度为，迭代精度为10-6。在命令中调用函数文件在命令中调用函数文件gauseidel.m，命令如下：，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=gauseidel(A,b,0,0,0,1.0e-6)8/2/2024电子信息学院例例
19、分别用分别用Jacobi迭代和迭代和Gauss-Serdel迭代法求解下列线性迭代法求解下列线性方程组，看是否收敛。方程组，看是否收敛。命令如下：命令如下：a=1,2,-2;1,1,1;2,2,1;b=9;7;6;x,n=jacobi(a,b,0;0;0)x,n=gauseidel(a,b,0;0;0)8/2/2024电子信息学院3.1.6 矩阵的秩与迹矩阵的秩与迹1矩阵的秩矩阵的秩矩阵线性无关的行数与列数称为矩阵的矩阵线性无关的行数与列数称为矩阵的秩。在秩。在MATLAB中，求矩阵秩的函数是中，求矩阵秩的函数是rank(A)。2矩阵的迹矩阵的迹矩阵的迹等于矩阵的对角线元素之和，矩阵的迹等于矩
20、阵的对角线元素之和，也等于矩阵的特征值之和。在也等于矩阵的特征值之和。在MATLAB中，求矩阵的迹的函数是中，求矩阵的迹的函数是trace(A)。8/2/2024电子信息学院3.1.7 向量和矩阵的范数向量和矩阵的范数矩阵或向量的范数用来度量矩阵或向量矩阵或向量的范数用来度量矩阵或向量在某种意义下的长度。范数有多种方法在某种意义下的长度。范数有多种方法定义，其定义不同，范数值也就不同。定义，其定义不同，范数值也就不同。8/2/2024电子信息学院1向量的向量的3种常用范数及其计算函数种常用范数及其计算函数在在MATLAB中，求向量范数的函数为：中，求向量范数的函数为：(1) norm(V)或或
21、norm(V,2)：计算向量：计算向量V的的2范数。范数。(2) norm(V,1)：计算向量：计算向量V的的1范数。范数。(3) norm(V,inf)：计算向量：计算向量V的的范数。范数。2矩阵的范数及其计算函数矩阵的范数及其计算函数MATLAB提供了求提供了求3种矩阵范数的函数，种矩阵范数的函数，其函数调用格式与求向量的范数的函数其函数调用格式与求向量的范数的函数完全相同。完全相同。8/2/2024电子信息学院3.1.8 矩阵的条件数矩阵的条件数在在MATLAB中，计算矩阵中，计算矩阵A的的3种条件数的种条件数的函数是：函数是：(1) cond(A,1) 计算计算A的的1范数下的条件数。
22、范数下的条件数。(2) cond(A)或或cond(A,2) 计算计算A的的2范数范数数下的条件数。数下的条件数。(3) cond(A,inf) 计算计算A的的 范数下的条范数下的条件数。件数。8/2/2024电子信息学院3.1.9 矩阵的特征值与特征向量矩阵的特征值与特征向量在在MATLAB中，计算矩阵中，计算矩阵A的特征值和的特征值和特征向量的函数是特征向量的函数是eig(A)，常用的调用格，常用的调用格式有式有2种：种：(1) E=eig(A)：求矩阵：求矩阵A的全部特征值，的全部特征值，构成向量构成向量E。(2) V,D=eig(A)：求矩阵：求矩阵A的全部特征的全部特征值，构成对角阵
23、值，构成对角阵D，并求，并求A的特征向量构的特征向量构成成V的列向量。的列向量。8/2/2024电子信息学院例例 用求特征值的方法解方程。用求特征值的方法解方程。3x5-7x4+5x2+2x-18=0p=3,-7,0,5,2,-18;A=compan(p); %A的伴随矩阵的伴随矩阵x1=eig(A) %求求A的特征值的特征值x2=roots(p) %直接求多项式直接求多项式p的零点的零点作业作业:p120第第8、14题题8/2/2024电子信息学院3.2 多项式计算多项式计算3.2.1 多项式的四则运算多项式的四则运算1多项式的加减运算（详见课本多项式的加减运算（详见课本p87）作业作业:编
24、写子函数可对任意二个多项式进行加减操作编写子函数可对任意二个多项式进行加减操作.(自动补零自动补零)2多项式乘法运算多项式乘法运算函数函数conv(P1,P2)用于求多项式用于求多项式P1和和P2的乘积。这里，的乘积。这里，P1、P2是两个多项式系数向量。是两个多项式系数向量。作业：求多项式作业：求多项式x4+8x3-10与多项式与多项式2x2-x+3的乘积。的乘积。提高提高:对多项式进行四则运算对多项式进行四则运算,输入表达式而不是输入表达式而不是向量向量,能输出运算结果能输出运算结果.8/2/2024电子信息学院3多项式除法多项式除法函数函数Q,r=deconv(P1,P2)用于对多项式用
25、于对多项式P1和和P2作除法运算。其中作除法运算。其中Q返回多项式返回多项式P1除以除以P2的商的商式，式，r返回返回P1除以除以P2的余式。这里，的余式。这里，Q和和r仍是仍是多项式系数向量。多项式系数向量。deconv是是conv的逆函数，即有的逆函数，即有P1=conv(P2,Q)+r。作业：作业： 求多项式求多项式x4+8x3-10除以多项式除以多项式2x2-x+3的的结果。结果。8/2/2024电子信息学院3.2.2 多项式的导函数多项式的导函数对多项式求导数的函数是：对多项式求导数的函数是：p=polyder(P)：求多项式：求多项式P的导函数的导函数p=polyder(P,Q)：
26、求：求PQ的导函数的导函数p,q=polyder(P,Q)：求：求P/Q的导函数，导的导函数，导函数的分子存入函数的分子存入p，分母存入，分母存入q。上述函数中，参数上述函数中，参数P,Q是多项式的向量表示，是多项式的向量表示，结果结果p,q也是多项式的向量表示。也是多项式的向量表示。8/2/2024电子信息学院例例 求有理分式的导数。求有理分式的导数。命令如下：命令如下：P=1;Q=1,0,5;p,q=polyder(P,Q)8/2/2024电子信息学院3.2.3 多项式的求值多项式的求值MATLAB提供了两种求多项式值的函数：提供了两种求多项式值的函数：polyval与与polyvalm，
27、它们的输入参数均为多项，它们的输入参数均为多项式系数向量式系数向量P和自变量和自变量x。两者的区别在于前者。两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。是代数多项式求值，而后者是矩阵多项式求值。8/2/2024电子信息学院1代数多项式求值代数多项式求值polyval函数用来求代数多项式的值，其调函数用来求代数多项式的值，其调用格式为：用格式为：Y=polyval(P,x)若若x为一数值，则求多项式在该点的值；若为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。个元素求其多项式的值。作业：已知多项式作业：
28、已知多项式x4+8x3-10，分别取，分别取x=1.2和一个和一个23矩阵为自变量计算该多矩阵为自变量计算该多项式的值。项式的值。8/2/2024电子信息学院2矩阵多项式求值矩阵多项式求值polyvalm函数用来求矩阵多项式的值，其调用格函数用来求矩阵多项式的值，其调用格式与式与polyval相同，但含义不同。相同，但含义不同。polyvalm函数函数要求要求x为方阵，它以方阵为自变量求多项式的为方阵，它以方阵为自变量求多项式的值。设值。设A为方阵，为方阵，P代表多项式代表多项式x3-5x2+8，那么，那么polyvalm(P,A)的含义是：的含义是：A*A*A-5*A*A+8*eye(siz
29、e(A)而而polyval(P,A)的含义是：的含义是：A.*A.*A-5*A.*A+8*ones(size(A)作业：作业： 仍以多项式仍以多项式x4+8x3-10为例，取一个为例，取一个22矩矩阵为自变量分别用阵为自变量分别用polyval和和polyvalm计算该多计算该多项式的值。项式的值。8/2/2024电子信息学院3.2.4 多项式求根多项式求根n次多项式具有次多项式具有n个根，当然这些根可能是个根，当然这些根可能是实根，也可能含有若干对共轭复根。实根，也可能含有若干对共轭复根。MATLAB提供的提供的roots函数用于求多项式函数用于求多项式的全部根，其调用格式为：的全部根，其调
30、用格式为：x=roots(P)其中其中P为多项式的系数向量，求得的根赋给为多项式的系数向量，求得的根赋给向量向量x，即，即x(1),x(2),x(n)分别代表多项分别代表多项式的式的n个根。个根。8/2/2024电子信息学院例例 求多项式求多项式x4+8x3-10的根。的根。命令如下：命令如下：A=1,8,0,0,-10;x=roots(A)若已知多项式的全部根，则可以用若已知多项式的全部根，则可以用poly函函数建立起该多项式，其调用格式为：数建立起该多项式，其调用格式为：P=poly(x)若若x为具有为具有n个元素的向量，则个元素的向量，则poly(x)建立建立以以x为其根的多项式，且将该
31、多项式的系为其根的多项式，且将该多项式的系数赋给向量数赋给向量P。8/2/2024电子信息学院例例 已知已知 f(x)(1) 计算计算f(x)=0 的全部根。的全部根。(2) 由方程由方程f(x)=0的根构造一个多项式的根构造一个多项式g(x)，并与，并与f(x)进行对比。进行对比。命令如下：命令如下：P=3,0,4,-5,-7.2,5;X=roots(P) %求方程求方程f(x)=0的根的根G=poly(X) %求多项式求多项式g(x)作业作业:p120第第4、5题题8/2/2024电子信息学院8/2/2024电子信息学院3.2.5 数据插值数据插值3.2.5.1 一维数据插值一维数据插值在
32、在MATLAB中，实现这些插值的函数是中，实现这些插值的函数是interp1，其调用格式为：，其调用格式为：Y1=interp1(X,Y,X1,method)函数根据函数根据X,Y的值，计算函数在的值，计算函数在X1处的值。处的值。X,Y是两个等长的已知向量，分别描述采是两个等长的已知向量，分别描述采样点和样本值，样点和样本值，X1是一个向量或标量，是一个向量或标量，描述欲插值的点，描述欲插值的点，Y1是一个与是一个与X1等长的等长的插值结果。插值结果。method是插值方法，允许的是插值方法，允许的取值有取值有linear、nearest、cubic、spline。8/2/2024电子信息学
33、院注意：注意：X1的取值范围不能超出的取值范围不能超出X的给定范的给定范围，否则，会给出围，否则，会给出“NaN”错误。错误。例用不同的插值方法计算在例用不同的插值方法计算在/2点的值。点的值。MATLAB中有一个专门的中有一个专门的3次样条插值函次样条插值函数数Y1=spline(X,Y,X1)，其功能及使用方，其功能及使用方法与函数法与函数Y1=interp1(X,Y,X1,spline)完完全相同。全相同。8/2/2024电子信息学院例某观测站测得某日例某观测站测得某日6:00时至时至18:00时之间每隔时之间每隔2小时的室小时的室内外温度内外温度()，用，用3次样条插值分别求得该日室内
34、外次样条插值分别求得该日室内外6:30至至17:30时之间每隔时之间每隔2小时各点的近似温度小时各点的近似温度()。设时间变量设时间变量h为一行向量，温度变量为一行向量，温度变量t为一个两列矩阵，其为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令中第一列存放室内温度，第二列储存室外温度。命令如下：如下：h =6:2:18;t=18,20,22,25,30,28,24;15,19,24,28,34,32,30;XI =6.5:2:17.5YI=interp1(h,t,XI,spline) %用用3次样条插值计算次样条插值计算8/2/2024电子信息学院3.2.5.2 二维数据插值
35、二维数据插值在在MATLAB中，提供了解决二维插值问题的函数中，提供了解决二维插值问题的函数interp2，其调用格式为：，其调用格式为：Z1=interp2(X,Y,Z,X1,Y1,method)其中其中X,Y是两个向量，分别描述两个参数的采样是两个向量，分别描述两个参数的采样点，点，Z是与参数采样点对应的函数值，是与参数采样点对应的函数值，X1,Y1是是两个向量或标量，描述欲插值的点。两个向量或标量，描述欲插值的点。Z1是根据是根据相应的插值方法得到的插值结果。相应的插值方法得到的插值结果。 method的的取值与一维插值函数相同。取值与一维插值函数相同。X,Y,Z也可以是矩也可以是矩阵形
36、式。阵形式。同样，同样，X1,Y1的取值范围不能超出的取值范围不能超出X,Y的给定范的给定范围，否则，会给出围，否则，会给出“NaN”错误。错误。8/2/2024电子信息学院例设例设z=x2+y2，对，对z函数在函数在0,10,2区域内进行插值。区域内进行插值。例某实验对一根长例某实验对一根长10米的钢轨进行热源的温度传播测试。米的钢轨进行热源的温度传播测试。用用x表示测量点表示测量点0:2.5:10(米米)，用，用h表示测量时间表示测量时间0:30:60(秒秒)，用，用T表示测试所得各点的温度表示测试所得各点的温度()。试用。试用线性插值求出在一分钟内每隔线性插值求出在一分钟内每隔20秒、钢
37、轨每隔秒、钢轨每隔1米处的米处的温度温度TI。命令如下：命令如下：x=0:2.5:10;h=0:30:60;T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41;xi=0:10;hi=0:20:60;TI=interp2(x,h,T,xi,hi)8/2/2024电子信息学院3.2.5.3 曲线拟合曲线拟合在在MATLAB中，用中，用polyfit函数来求得最小二乘拟函数来求得最小二乘拟合多项式的系数，再用合多项式的系数，再用polyval函数按所得的多函数按所得的多项式计算所给出的点上的函数近似值。项式计算所给出的点上的函数近似值。polyfit函数的调用格式为
38、：函数的调用格式为：P,S=polyfit(X,Y,m)函数根据采样点函数根据采样点X和采样点函数值和采样点函数值Y，产生一个，产生一个m次多项式次多项式P及其在采样点的误差向量及其在采样点的误差向量S。其中。其中X,Y是两个等长的向量，是两个等长的向量，P是一个长度为是一个长度为m+1的的向量，向量，P的元素为多项式系数。的元素为多项式系数。polyval函数的功能是按多项式的系数计算函数的功能是按多项式的系数计算x点多点多项式的值，将在项式的值，将在6.5.3节中详细介绍。节中详细介绍。8/2/2024电子信息学院例已知数据表例已知数据表t,y，试求，试求2次拟合多项式次拟合多项式p(t)
39、，然后求，然后求ti=1,1.5,2,2.5,9.5,10各点的函数近似值。各点的函数近似值。8/2/2024电子信息学院3.3 数据统计处理数据统计处理3.3.1 最大值和最小值最大值和最小值MATLAB提供的求数据序列的最大值和最提供的求数据序列的最大值和最小值的函数分别为小值的函数分别为max和和min，两个函数，两个函数的调用格式和操作过程类似。的调用格式和操作过程类似。1求向量的最大值和最小值求向量的最大值和最小值求一个向量求一个向量X的最大值的函数有两种调用的最大值的函数有两种调用格式，分别是：格式，分别是：(1) y=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，
40、如果，如果X中包含复数元素，则按模取最大中包含复数元素，则按模取最大值。值。8/2/2024电子信息学院(2) y,I=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，最大值的序号存入最大值的序号存入I，如果，如果X中包含复数元素，中包含复数元素，则按模取最大值。则按模取最大值。求向量求向量X的最小值的函数是的最小值的函数是min(X)，用法和，用法和max(X)完全相同。完全相同。例例6-1 求向量求向量x的最大值。的最大值。命令如下：命令如下：x=-43,72,9,16,23,47;y=max(x) %求向量求向量x中的最大值中的最大值y,l=max(x) %求向量求向量x中
41、的最大值及其该元中的最大值及其该元素的位置素的位置8/2/2024电子信息学院2求矩阵的最大值和最小值求矩阵的最大值和最小值求矩阵求矩阵A的最大值的函数有的最大值的函数有3种调用格式，种调用格式，分别是：分别是：(1) max(A)：返回一个行向量，向量的第：返回一个行向量，向量的第i个元素是矩阵个元素是矩阵A的第的第i列上的最大值。列上的最大值。(2) Y,U=max(A)：返回行向量：返回行向量Y和和U，Y向量记录向量记录A的每列的最大值，的每列的最大值，U向量记录向量记录每列最大值的行号。每列最大值的行号。8/2/2024电子信息学院(3) max(A,dim)：dim取取1或或2。di
42、m取取1时，时，该函数和该函数和max(A)完全相同；完全相同；dim取取2时，时，该函数返回一个列向量，其第该函数返回一个列向量，其第i个元素是个元素是A矩阵的第矩阵的第i行上的最大值。行上的最大值。求最小值的函数是求最小值的函数是min，其用法和，其用法和max完全完全相同。相同。例例6-2 分别求分别求34矩阵矩阵x中各列和各行元素中各列和各行元素中的最大值，并求整个矩阵的最大值和中的最大值，并求整个矩阵的最大值和最小值。最小值。8/2/2024电子信息学院3两个向量或矩阵对应元素的比较两个向量或矩阵对应元素的比较函数函数max和和min还能对两个同型的向量或矩阵进还能对两个同型的向量或
43、矩阵进行比较，调用格式为：行比较，调用格式为：(1) U=max(A,B)：A,B是两个同型的向量或矩阵，是两个同型的向量或矩阵，结果结果U是与是与A,B同型的向量或矩阵，同型的向量或矩阵，U的每个元的每个元素等于素等于A,B对应元素的较大者。对应元素的较大者。(2) U=max(A,n)：n是一个标量，结果是一个标量，结果U是与是与A同同型的向量或矩阵，型的向量或矩阵，U的每个元素等于的每个元素等于A对应元对应元素和素和n中的较大者。中的较大者。min函数的用法和函数的用法和max完全相同。完全相同。例例6-3 求两个求两个23矩阵矩阵x, y所有同一位置上的较大所有同一位置上的较大元素构成
44、的新矩阵元素构成的新矩阵p。8/2/2024电子信息学院3.3.2 求和与求积求和与求积数据序列求和与求积的函数是数据序列求和与求积的函数是sum和和prod，其使用方法类似。设，其使用方法类似。设X是一个向量，是一个向量，A是一个矩阵，函数的调用格式为：是一个矩阵，函数的调用格式为：sum(X)：返回向量：返回向量X各元素的和。各元素的和。prod(X)：返回向量：返回向量X各元素的乘积。各元素的乘积。sum(A)：返回一个行向量，其第：返回一个行向量，其第i个元素是个元素是A的第的第i列的元素和。列的元素和。8/2/2024电子信息学院prod(A)：返回一个行向量，其第：返回一个行向量，
45、其第i个元素个元素是是A的第的第i列的元素乘积。列的元素乘积。sum(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于sum(A)；当；当dim为为2时，返回一个列向量，时，返回一个列向量，其第其第i个元素是个元素是A的第的第i行的各元素之和。行的各元素之和。prod(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于prod(A)；当；当dim为为2时，返回一个列向量，时，返回一个列向量，其第其第i个元素是个元素是A的第的第i行的各元素乘积。行的各元素乘积。例例6-4 求矩阵求矩阵A的每行元素的乘积和全部元的每行元素的乘积和全部元素的乘积。素的乘积。8/2/20
46、24电子信息学院3.3.3 平均值和中值平均值和中值求数据序列平均值的函数是求数据序列平均值的函数是mean，求数据序列中值的函数是，求数据序列中值的函数是median。两。两个函数的调用格式为：个函数的调用格式为：mean(X)：返回向量：返回向量X的算术平均值。的算术平均值。median(X)：返回向量：返回向量X的中值。的中值。mean(A)：返回一个行向量，其第返回一个行向量，其第i个元素是个元素是A的第的第i列的算术平均值。列的算术平均值。median(A)：返回一个行向量，其第返回一个行向量，其第i个元素是个元素是A的第的第i列的中值。列的中值。mean(A,dim)：当当dim为
47、为1时，该函数等同于时，该函数等同于mean(A)；当；当dim为为2时，返回一个列向量，时，返回一个列向量，其第其第i个元素是个元素是A的第的第i行的算术平均值。行的算术平均值。median(A,dim)：当当dim为为1时，该函数等同于时，该函数等同于median(A)；当；当dim为为2时，返回一个列向量，时，返回一个列向量，其第其第i个元素是个元素是A的第的第i行的中值。行的中值。例例6-5 分别求向量分别求向量x与与y的平均值和中值。的平均值和中值。8/2/2024电子信息学院3.3.4 累加和与累乘积累加和与累乘积在在MATLAB中，使用中，使用cumsum和和cumprod函数能
48、方便地求得向量和函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的调用格式为：矩阵元素的累加和与累乘积向量，函数的调用格式为：cumsum(X)：返回向量：返回向量X累加和向量。累加和向量。cumprod(X)：返回向量：返回向量X累乘积向量。累乘积向量。cumsum(A)：返回一个矩阵，其第：返回一个矩阵，其第i列是列是A的第的第i列的累加和向量。列的累加和向量。cumprod(A)：返回一个矩阵，其第：返回一个矩阵，其第i列是列是A的第的第i列的累乘积向量。列的累乘积向量。cumsum(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于cumsum(A)；当；当dim
49、为为2时，返回一个矩阵，其第时，返回一个矩阵，其第i行是行是A的第的第i行的累加和向量。行的累加和向量。cumprod(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于cumprod(A)；当；当dim为为2时，返回一个向量，其第时，返回一个向量，其第i行是行是A的第的第i行的累乘积向量。行的累乘积向量。例例6-6 求求s的值。的值。8/2/2024电子信息学院3.3.5 标准方差与相关系数标准方差与相关系数1求标准方差求标准方差在在MATLAB中，提供了计算数据序列的标准方差的函数中，提供了计算数据序列的标准方差的函数std。对于向量。对于向量X，std(X)返回一个标准方差
50、。对于矩阵返回一个标准方差。对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵返回一个行向量，它的各个元素便是矩阵A各列或各行的标准方差。各列或各行的标准方差。std函数的一般调用格式为：函数的一般调用格式为：Y=std(A,flag,dim)其中其中dim取取1或或2。当。当dim=1时，求各列元素的标准方差；时，求各列元素的标准方差；当当dim=2时，则求各行元素的标准方差。时，则求各行元素的标准方差。flag取取0或或1，当当flag=0时，按时，按1所列公式计算标准方差，当所列公式计算标准方差，当flag=1时，时，按按2所列公式计算标准方差。缺省所列公式计算标准方差。缺省f
51、lag=0，dim=1。例例6-7 对二维矩阵对二维矩阵x，从不同维方向求出其标准方差。，从不同维方向求出其标准方差。8/2/2024电子信息学院2相关系数相关系数MATLAB提供了提供了corrcoef函数，可以求出函数，可以求出数据的相关系数矩阵。数据的相关系数矩阵。corrcoef函数的调函数的调用格式为：用格式为：corrcoef(X)：返回从矩阵：返回从矩阵X形成的一个相关形成的一个相关系数矩阵。此相关系数矩阵的大小与矩系数矩阵。此相关系数矩阵的大小与矩阵阵X一样。它把矩阵一样。它把矩阵X的每列作为一个变的每列作为一个变量，然后求它们的相关系数。量，然后求它们的相关系数。corrco
52、ef(X,Y)：在这里，：在这里，X,Y是向量，它们是向量，它们与与corrcoef(X,Y)的作用一样。的作用一样。8/2/2024电子信息学院例例6-8 生成满足正态分布的生成满足正态分布的100005随机矩随机矩阵，然后求各列元素的均值和标准方差，阵，然后求各列元素的均值和标准方差，再求这再求这5列随机数据的相关系数矩阵。列随机数据的相关系数矩阵。命令如下：命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)8/2/2024电子信息学院3.3.6 排序排序MATLAB中对向量中对向量X是排序函数是是排序函数是sort(X)，函数，函数返
53、回一个对返回一个对X中的元素按升序排列的新向量。中的元素按升序排列的新向量。sort函数也可以对矩阵函数也可以对矩阵A的各列或各行重新排序，的各列或各行重新排序，其调用格式为：其调用格式为：Y,I=sort(A,dim)其中其中dim指明对指明对A的列还是行进行排序。若的列还是行进行排序。若dim=1，则按列排；若，则按列排；若dim=2，则按行排。，则按行排。Y是排序是排序后的矩阵，而后的矩阵，而I记录记录Y中的元素在中的元素在A中位置。中位置。例例6-9 对二维矩阵做各种排序。对二维矩阵做各种排序。8/2/2024电子信息学院3.4 函数分析与数值积分函数分析与数值积分3.4.1 函数在函
54、数在MATLAB中的表示与函数的绘中的表示与函数的绘1.函数表示与计算函数表示与计算a.函数文件函数文件 p105页的页的humps.m表示表示b.Inline内联函数实现内联函数实现2.函数的绘制函数的绘制a.单变量函数绘制图形命令单变量函数绘制图形命令fplot函数的调用格式为：函数的调用格式为： fplot(fname,lims,tol,选项选项)其中其中fname为函数名，以字符串形式出现，为函数名，以字符串形式出现，lims为为x,y的取值范围，的取值范围，tol为相对允许误差，其系统默认值为为相对允许误差，其系统默认值为2e-3。选项定义与。选项定义与plot函数相同。函数相同。例
55、例5-9 用用fplot函数绘制函数绘制f(x)=cos(tan(x)的曲线。的曲线。命令如下：命令如下：fplot(cos(tan(pi*x), 0,1,1e-4)8/2/2024电子信息学院b.简易的函数绘图命令简易的函数绘图命令ezplotMATLAB提供了一个提供了一个ezplot函数绘制函数图形，函数绘制函数图形，下面介绍其用法。下面介绍其用法。(1) 对于函数对于函数f = f(x)，ezplot函数的调用格式为：函数的调用格式为：ezplot(f)：在默认区间：在默认区间-2x2绘制绘制f = f(x)的图的图形。形。ezplot(f, a,b)：在区间：在区间axb绘制绘制f
56、= f(x)的图形。的图形。8/2/2024电子信息学院(2) 对于函数对于函数f = f(x,y)，ezplot函数的调用格式为：函数的调用格式为：ezplot(f)：在默认区间：在默认区间-2x2和和-2y2绘制绘制f(x,y) = 0的图形。的图形。ezplot(f, xmin,xmax,ymin,ymax)：在区间：在区间xminxxmax和和yminyymax绘制绘制f(x,y) = 0的的图形。图形。ezplot(f, a,b)：在区间：在区间axb和和ay b绘制绘制f(x,y) = 0的图形。的图形。8/2/2024电子信息学院(3) 对于参数方程对于参数方程x = x(t)和
57、和y = y(t)，ezplot函数函数的调用格式为：的调用格式为：ezplot(x,y)：在默认区间：在默认区间0t2绘制绘制x=x(t)和和y=y(t)的图形。的图形。ezplot(x,y, tmin,tmax)：在区间：在区间tmin t tmax绘制绘制x=x(t)和和y=y(t)的图形。的图形。8/2/2024电子信息学院3.4.2 函数极值函数极值 MATLAB提供了基于单纯形算法求解函数极值提供了基于单纯形算法求解函数极值的函数的函数fmin和和fmins， 它们分别用于单变量函它们分别用于单变量函数和多变量函数的最小值，其调用格式为：数和多变量函数的最小值，其调用格式为： x=
58、fmin(fname,x1,x2) x=fmins(fname,x0)这两个函数的调用格式相似。其中这两个函数的调用格式相似。其中fmin函数用于函数用于求单变量函数的最小值点。求单变量函数的最小值点。fname是被最小化是被最小化的目标函数名，的目标函数名，x1和和x2限定自变量的取值范围。限定自变量的取值范围。fmins函数用于求多变量函数的最小值点，函数用于求多变量函数的最小值点，x0是求解的初始值向量。是求解的初始值向量。8/2/2024电子信息学院MATLAB没有专门提供求函数最大值的函数，但只要注没有专门提供求函数最大值的函数，但只要注意到意到-f(x)在区间在区间(a,b)上的最
59、小值就是上的最小值就是f(x)在在(a,b)的最的最大值，所以大值，所以fmin(f,x1,x2)返回函数返回函数f(x)在区间在区间(x1,x2)上上的最大值。的最大值。 例例 求求f(x)=x3-2x-5在在0,5内的最小值点。内的最小值点。 (1) 建立函数文件建立函数文件mymin.m。function fx=mymin(x)fx=x.3-2*x-5; (2) 调用调用fmin函数求最小值点。函数求最小值点。x=fmin(mymin,0,5)x= 0.81658/2/2024电子信息学院3.4.2数值积分与微分数值积分与微分3.4.2.1 数值积分基本原理数值积分基本原理 求解定积分的
60、数值方法多种多样，如简求解定积分的数值方法多种多样，如简单的梯形法、辛普生单的梯形法、辛普生(Simpson) 法、牛法、牛顿柯特斯顿柯特斯(Newton-Cotes)法等都是经常法等都是经常采用的方法。它们的基本思想都是将整采用的方法。它们的基本思想都是将整个积分区间个积分区间a,b分成分成n个子区间个子区间xi,xi+1，i=1,2,n，其中，其中x1=a，xn+1=b。这样求。这样求定积分问题就分解为求和问题。定积分问题就分解为求和问题。8/2/2024电子信息学院3.4.2.2 数值积分的实现方法数值积分的实现方法1变步长辛普生法变步长辛普生法基于变步长辛普生法，基于变步长辛普生法，M
61、ATLAB给出了给出了quad函函数来求定积分。该函数的调用格式为：数来求定积分。该函数的调用格式为： I,n=quad(fname,a,b,tol,trace)其中其中fname是被积函数名。是被积函数名。a和和b分别是定积分的分别是定积分的下限和上限。下限和上限。tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取控制是否展现积分过程，若取非非0则展现积分过程，取则展现积分过程，取0则不展现，缺省时取则不展现，缺省时取trace=0。返回参数。返回参数I即定积分值，即定积分值，n为被积函数为被积函数的调用次数。的调用次数。8/2
62、/2024电子信息学院 例求定积分。例求定积分。 (1) 建立被积函数文件建立被积函数文件fesin.m。function f=fesin(x)f=exp(-0.5*x).*sin(x+pi/6); (2) 调用数值积分函数调用数值积分函数quad求定积分。求定积分。S,n=quad(fesin,0,3*pi)S = 0.9008n = 778/2/2024电子信息学院2牛顿柯特斯法牛顿柯特斯法基于牛顿柯特斯法，基于牛顿柯特斯法，MATLAB给出了给出了quad8函数来求定积分。该函数的调用格函数来求定积分。该函数的调用格式为：式为：I,n=quad8(fname,a,b,tol,trace)
63、其中参数的含义和其中参数的含义和quad函数相似，只是函数相似，只是tol的缺省值取的缺省值取10-6。 该函数可以更精确地该函数可以更精确地求出定积分的值，且一般情况下函数调求出定积分的值，且一般情况下函数调用的步数明显小于用的步数明显小于quad函数，从而保证函数，从而保证能以更高的效率求出所需的定积分值。能以更高的效率求出所需的定积分值。8/2/2024电子信息学院例求定积分。例求定积分。(1) 被积函数文件被积函数文件fx.m。function f=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);(2) 调用函数调用函数quad8求定积分。求定积分。I=quad
64、8(fx,0,pi)I = 2.46748/2/2024电子信息学院例分别用例分别用quad函数和函数和quad8函数求定积分的近似值，函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。并在相同的积分精度下，比较函数的调用次数。调用函数调用函数quad求定积分：求定积分：format long;fx=inline(exp(-x);I,n=quad(fx,1,2.5,1e-10)I = 0.28579444254766n = 658/2/2024电子信息学院 调用函数调用函数quad8求定积分：求定积分：format long;fx=inline(exp(-x);I,n=quad8
65、(fx,1,2.5,1e-10)I = 0.28579444254754n = 338/2/2024电子信息学院3被积函数由一个表格定义被积函数由一个表格定义在在MATLAB中，对由表格形式定义的函数关系的中，对由表格形式定义的函数关系的求定积分问题用求定积分问题用trapz(X,Y)函数。其中向量函数。其中向量X,Y定义函数关系定义函数关系Y=f(X)。例例 用用trapz函数计算定积分。函数计算定积分。命令如下：命令如下：X=1:0.01:2.5;Y=exp(-X); %生成函数关系数据向量生成函数关系数据向量trapz(X,Y)ans = 0.285796824163938/2/2024
66、电子信息学院3.4.2 二重定积分的数值求解二重定积分的数值求解3.4.2使用使用MATLAB提供的提供的dblquad函数函数就可以直接求出上述二重定积分的就可以直接求出上述二重定积分的数值解。该函数的调用格式为：数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求该函数求f(x,y)在在a,bc,d区域上区域上的二重定积分。参数的二重定积分。参数tol，trace的用的用法与函数法与函数quad完全相同。完全相同。8/2/2024电子信息学院例计算二重定积分例计算二重定积分(1) 建立一个函数文件建立一个函数文件fxy.m：function f=
67、fxy(x,y)global ki;ki=ki+1; %ki用于统计被积函数的调用次数用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2) 调用调用dblquad函数求解。函数求解。global ki;ki=0;I=dblquad(fxy,-2,2,-1,1)kiI = 1.57449318974494ki = 10388/2/2024电子信息学院3.4.3.2数值微分数值微分1 数值差分与差商数值差分与差商2 数值微分的实现数值微分的实现在在MATLAB中，没有直接提供求数值导数的函数，只有计中，没有直接提供求数值导数的函数，只有计算向前差分的函数算向前差分的
68、函数diff，其调用格式为：，其调用格式为：DX=diff(X)：计算向量：计算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n)：计算：计算X的的n阶向前差分。例如，阶向前差分。例如，diff(X,2)=diff(diff(X)。DX=diff(A,n,dim)：计算矩阵：计算矩阵A的的n阶差分，阶差分，dim=1时时(缺省缺省状态状态)，按列计算差分；，按列计算差分；dim=2，按行计算差分。，按行计算差分。8/2/2024电子信息学院例例 生成以向量生成以向量V=1,2,3,4,5,6为基础的范得蒙矩为基础的范得蒙矩阵，按列进行
69、差分运算。阵，按列进行差分运算。命令如下：命令如下：V=vander(1:6)DV=diff(V) %计算计算V的一阶差分的一阶差分8/2/2024电子信息学院例例 用不同的方法求函数用不同的方法求函数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,
70、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,.,x,gx,-); %作图作图8/2/2024电子信息学院3.4.3 常微分方程初值问题的数值解法常微分方程初值问题的数值解法3.4.3.1 龙格库塔法简介龙格库塔法简介3.4.3.
71、2 龙格库塔法的实现龙格库塔法的实现 基于龙格库塔法，基于龙格库塔法，MATLAB提供了求常微分方程数提供了求常微分方程数值解的函数，一般调用格式为：值解的函数，一般调用格式为： t,y=ode23(fname,tspan,y0) t,y=ode45(fname,tspan,y0)其中其中fname是定义是定义f(t,y)的函数文件名，该函数文件必须返的函数文件名，该函数文件必须返回一个列向量。回一个列向量。tspan形式为形式为t0,tf,表示求解区间。表示求解区间。y0是初始状态列向量。是初始状态列向量。t和和y分别给出时间向量和相应分别给出时间向量和相应的状态向量。的状态向量。8/2/2
72、024电子信息学院例设有初值问题，试求其数值解，并与精确解相比较例设有初值问题，试求其数值解，并与精确解相比较(精精确解为确解为y(t)=)。 (1) 建立函数文件建立函数文件funt.m。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2) 求解微分方程。求解微分方程。t0=0;tf=10;y0=2;t,y=ode23(funt,t0,tf,y0); %求数值解求数值解y1=sqrt(t+1)+1; %求精确解求精确解tyy1 y为数值解，为数值解，y1为精确值，显然两者近似。为精确值，显然两者近似。8/2/2024电子信息学院例例 求解著名的求解著名的Van der Pol方程。方程。例有例有Lorenz模型的状态方程，试绘制系统相平面图。模型的状态方程，试绘制系统相平面图。8/2/2024电子信息学院
