第3章 MATLAB在高等数学中的应用 3.1 矩阵分析 3.2 多项式运算 3.3 数据的分析与统计 3.4 函数分析与数值积分

3.1 矩阵分析 1．矢量范数和矩阵范数 矩阵范数是对矩阵的一种测度。矢量的p范数和矩阵A的p范数分别定为： 当p＝2时为常用的欧拉范数，一般p还可取l和∞。这在MATLAB中可利用norm函数实现，p缺省时为p=2。 格式：n＝norm(A) 功能：计算矩阵A的最大奇异值，相当于n=max(svd(A))。 格式：n＝norm(A，p) 功能：norm函数可计算几种不同类型的矩阵范数，根据p的不同可得到不同的范数

⑴矩阵求逆函数inv及行列式值函数det 2．矩阵求逆及行列式值 ⑴矩阵求逆函数inv及行列式值函数det 逆矩阵的定义：对于任意阶 n×n 方阵A，如果能找到一个同阶的方阵V，使得满足：A*V=I。其中I为n阶的单位矩阵eye(n)。则V就是A的逆矩阵。数学符号表示为：V=A-1。逆矩阵V存在的条件是A的行列式不等于0。 格式：V=inv(A) 功能：返回方阵A的逆矩阵V。 格式：X=det(A) 功能：计算方阵A的行列式值。 ⑵伪逆矩阵函数pinv 伪逆矩阵的MATLAB定义：从数学意义上讲，当矩阵A为非方阵时，其矩阵的逆是不存在的。在MATLAB中，为了求线性方程组的需要，把inv(A′*A)*A′的运算定义为伪逆函数pinv，这样对非方阵，利用伪逆函数pinv可以求得矩阵的伪逆，伪逆在一定程度上代表着矩阵的逆。 格式：C=pinv(A) 功能：计算非方阵A的伪逆矩阵。

3．线性代数方程求解 一般线性方程组的 写成矩阵形式可表示为：AX＝B 或 XA＝B。其中系数矩阵A的阶数为m×n。在MATLAB中，引入矩阵除法求解。 (1)求解方程AX=B 格式：X=A\B 条件：矩阵A与矩阵B的行数必须相等。 (2)求解方程XA=B 格式： X=B/A 条件：矩阵A与矩阵B的列数必须相等。

4．矩阵的分解 (1)三角(LU)分解函数lu 所谓三角解就是将一个方阵表示成两个基本三角阵的乘积(A=LU），其中一个为下三角矩阵L，另一个为上三角形矩阵U，因而矩阵的三角分解又叫LU分解或叫LR分解。矩阵 分解的两个矩阵分别可表示为： 格式一：[L,U]=lu(A) 功能：返回一个上三角矩阵U和一个置换下三角矩阵L(即下三角矩阵与置换矩阵的乘积)，满足A=L*U。 格式二：[L,U,P]=lu(A) 功能：返回上三角矩阵U，真正下三角矩阵L，及一个置换矩阵P(用来表示排列规则的矩阵)，满足L*U=P*A；如果P为单位矩阵，满足A=L*U。

(2)正交(QR)分解函数 将矩阵A分解为一个正交矩阵与另一个矩阵的乘积称为矩阵A的正交分解。 格式一：[Q, R]=qr(A) 功能：产生与A同维的上三角矩阵R和一个实正交矩阵或复归一化矩阵Q，满足：A=Q*R，Q’*Q=I。 格式二：[Q,R,E]=qr(A) 功能：产生一个置换矩阵E，一个上三角矩阵R(其对角线元素降序排列)和一个归一化矩阵Q，满足A*E=Q*R；

矩阵A的奇异值和相应的一对奇异矢量u、v满足： 5．奇异值分解 矩阵A的奇异值和相应的一对奇异矢量u、v满足： 同样利用奇异值构成对角阵，相应的奇异矢量作为列构成两个正交矩阵U、V，则有： 其中AT表示转置矩阵。由于U和V正交，因此可得奇异值分解： 格式一：[U,S,V]=svd(x) 功能：返回3个矩阵，使得X=U*S*V’。其中S为与X相同维数的矩阵，且其对角元素为非负递减。 格式二：S=svd(A) 功能：返回奇异值组成的向量。

6．矩阵的特征值分析 矩阵A的特征值 和特征矢量 ，满足： 以特征值构成对角阵 ，相应的特征矢量作为列构成矩阵V，则有： 格式一：d=eig(A) 功能：返回方阵A的全部特征值所构成的向量。 格式二：[V,D]=eig(A) 功能：返回矩阵V和D。其中对角阵D的对角元素为A的特征值，V的列向量是相应的特征向量，使得A*V=V*D。

7．矩阵的幂次运算: A^p 在MATLAB中，矩阵的幂次运算是指以下两种情况： 1、矩阵为底数，指数是标量的运算操作； 2、底数是标量，矩阵为指数的运算操作。 两种情况都要求矩阵是方阵，否则，将显示出错信息。 (1) 矩阵的正整数幂 如果A是一个方阵，p是一个正整数，那么幂次表示A自己乘p次。 (2) 矩阵的负数幂 如果A是一个非奇异方阵，p是一个正整数，那么A^(－p)表示inv(A)自己乘p次。 (3) 矩阵的分数幂 如果A是一个方阵，p取分数，它的结果取决于矩阵的特征值的分布。 (4) 矩阵的元素幂、按矩阵元素的幂 利用运算符“A.^p”实现矩阵的元素幂或按矩阵元素的幂运算。

8．矩阵结构形式的提取与变换 (1) 矩阵左右翻转函数fliplr( ) 格式：X=fliplr(A) (2) 矩阵上下翻转函数flipud 格式：X=flipud(A) (3) 矩阵阶数重组函数reshape 格式一：X=reshape(A,n,m) 功能：将矩阵A中的所有元素按列的秩序重组成n×m阶矩阵X，当A中没有m×n个元素时会显示出错信息。 格式二 ：X=reshape(A,m,n,p,...) 或 X=reshape(A,[m,n,p,...]) 功能：从A中形成多维阵列(m×n×p×...)。

功能：将矩阵按反时针旋转k*90o，其中k应为整数。 (4) 矩阵整体反时针旋转函数rot90( ) 格式一：X=rot90(A) 功能：将矩阵按反时针旋转90o。 格式二：X=rot90(A, k) 功能：将矩阵按反时针旋转k*90o，其中k应为整数。 (5) 对角矩阵和矩阵的对角化函数diag( ) 格式一：X=diag(A,k) 功能：当A为n元向量时，可得n+abs(k)阶的方阵X，其A的元素处于第k条对角线上；k=0表示主对角线，k>0表示在主对角线之上，k<0表示在主对角线之下。当A为矩阵时，X=diag(A,k)得到列向量X，它取自于X的第k个对角线上的元素。 格式二：X=diag(A) 功能：当A为n元向量时，等同于k=0时的X=diag(A,k)，即产生A的元素处于主对角线的对角方阵。当A为矩阵时，X=diag(A)相当于k=0。

(6) 取矩阵的左下三角部分函数tril( ) 格式一：X=tril(A,k) 功能：得到矩阵A的第k条对角线及其以下的元素；当k=0时表示主对角线，k>0表示主对角线之上，k<0表示主对角线以下。 格式二：X=tril(A) 功能：得到矩阵A的下三角阵。 (7) 取矩阵的右上三角部分函数triu( ) 格式一：X=triu(A,k) 功能：得到矩阵A的第k条对角线及其以上的元素；当k=0时表示主对角线，k>0表示主对角线之上，k<0表示主对角线以下。 格式二：X=triu(A) 功能：得到矩阵A的右上三角阵。 (8) 利用“：”将矩阵元素按列取出排成一列 方法：X=A(:)’

3.2 多项式运算 3.2.1 多项式表示及其四则运算 1．MATLAB的多项式表示 对多项式： 用其系数的行向量p=[an, an-1, …,a1, a0]来表示。注意：如果 x的某次幂的系数为零，这个零必须列入系数向量中。例如一个一元3次多项式： 可表示成行向量：p=[1,0,–2, 5]。 2．多项式的加减运算 格式：A=B+C 3．多项式相乘运算 格式：w=conv(u,v) 功能：返回u和v两向量的卷积，也就是u和v代表的两多项式的乘积。 4．多项式相除 格式：[q , r]=deconv(u , v) 功能：给出商多项式q和余数多项式r ，u为被除多项式

3.2.2 多项式求导、求根和求值 1．多项式求导函数polyder 格式一：k=polyder(p) 功能：返回多项式p的一阶导数。 格式二：k=polyder(u,v) 功 能：返回多项式u与v乘积的导数。 格式三：[q,d]=polyer(u,v) 功 能：返回多项式商u/v的导数, 返回的格式为：q为分子, d为分母。 2．多项式的根 求解多项式的根，即p(x)=0的解。 格式：r=roots(p) 功能：返回多项式p(x)的根。注意，MATLAB按惯例规定，多项式是行向量，根是列向量。

3．多项式求值函数polyval( ) 利用函数polyval可以求得多项式在某一点的值。 格式：y=polyval(p,x) 功能：返回多项式p在x处的值。其中x可以是复数，也可以是数组。 当多项式的变量是矩阵时，构成的矩阵多项式可以利用polyvalm函数求值。 格式：Y=polyvalm(p,X) 功能：返回矩阵多项式p在X处的值。 4．部分分式展开函数residue( ) 格式一：[r,p,k]=residue(b,a) 功能：把b(s)/a(s)展开成： 其中r代表余数数组，p代表极点数组，ks代表部分分式展开的常数项。当分母多项式的阶次数高于分子多项式的阶次数时ks=0 格式二：[b, a]=residue(r, p, k) 功能：格式一的逆作用。

3.3.3 多项式拟合与多项式插值 1．多项式拟合函数polyfit( ) 格式：p=polyfit(x,y,n) 功能：利用已知的数据向量x和y所确定的数据点，采用最小二乘法构造出n阶多项式去逼近已知的离散数据，实现多项式曲线的拟合。其中p是求出的多项式系数，n阶多项式应该有n+1个系数，故p的长度为n+1。 2．多项式插值 插值和拟合的不同点在于：①插值函数通常是分段的，人们关心的不是函数的表达式，而是插值出的数据点；②插值函数应通过给定的数据点。 (1)一维插值函数interpl( ) 格式：yi＝interpl (x, y, xi, 'method') 功能：为给定的数据对(x,y)以及x坐标上的插值范围向量xi，用指定所使用的插值方法method实现插值。yi是插值后的对应数据点集的y坐标。插值的方法method有以下6种可供选择： nearest(最邻近插值法)、linear(线性插值)、 spline(三次样条插值)、cubic(立方插值)、pchip(三次Hermite插值)、v5cubic。

（2）二维插值函数 格式：zi=griddata(x,y,z,xi,yi,method) 功能：非等距插值。已知的元素值由3个向量来描述：x、y和z。函数返回值为一矩阵zi，其元素的值由x、y和z确定的二元函数插值得到。method可为：‘linear’(默认值)、‘cubic’ 、‘nearest’ 、‘v4’ 。 格式：zi=interp2(x, y, z, xi, yi, 'method') 功能：单调节点插值。已知的元素值由3个向量来描述：x、y和z。其中，x、y是已知数据组并且大小相同，z是相对应的已知点上的函数值；xi、yi是用于插值的矢量； zi是根据相应的插值方法并且与(xi, yi)对应的插值结果。插值方法method：‘linear’(默认值)、‘cubic’ 、‘nearest’ 、’spline’ 。

(3) 高维插值和交互式样条插值 ·高维插值函数：三维插值及三维以上的插值称为高维插值。 用于实现高维插值的函数有：interp3(三维插值函数)、interpn(n维插值函数)、ndgrid(n维数据网格)。其调用格式与interp2函数很类似，这里就不再重述了，具体调用格式读者可利用help来得。 ·交互式样条插值函数 在MATLAB 6.0及以上版本中，样条工具箱新增加了交互式插值样条函数splinetool。该函数以对话框的形式为用户提供了插值过程。 格式一：splinetool 功能：用于生成各种样条曲线，这里几乎包括所有生成样条曲线方法。在它的初始菜单中提供了各种数据，用户可以选择一种生成的样条曲线。 格式二：splinetool(x,y) 功能：用户输入数组x、y，并在用户图形界面下生成样条曲线。

3.3 数据分析与统计 3.3.1 数据基本操作 1．求最大值函数max 格式一：xM＝max(x) 3.3 数据分析与统计 3.3.1 数据基本操作 1．求最大值函数max 格式一：xM＝max(x) 功能：如果x是向量，返回x中最大值元素；如果x是矩阵，则将矩阵每列作为处理向量，返回一个行向量，其元素为矩阵每列中的最大元素；如果x为多维数组，则沿第一个非单元素维进行处理，求得各向量的最大值。 格式二：xM=max(x, y) 功能：返回一个与x和y一样大小的数组，其元素取x或y中最大的一个。 格式三：xM=max(x，[]，dim) 功能：返回数组(矩阵)x由标量dim所指定的维数(或行)中的最大值。 格式四：[xM，I]＝max(…) 功能：返回最大值同时，返回一个下标向量。 如果输入数据x为复数，max函数返回复数最大模：max(abs(x))。

2．求最小值函数min min函数的调用格式与max函数的调用格式相同，只是功能与max函数相反，所得结果为最小值。如果输入数据x为复数，min函数返回复数最小模：min(abs(x)) 3．求平均值函数mean 格式：M＝mean(x) 功能：如果x为向量，则返回向量x的平均值；如果x为矩阵，则将矩阵每列当作向量处理，返回一个平均值行向量；如果A为多维数组，则沿第一个非单元素维进行处理，返回一个平均值数组。 4. 求中间值函数median 格式：M=median(x) 功能：如果x为向量，则返回向量x的中间值；如果x为矩阵，则将矩阵每列当作为处理向量，返回一个中间值行向量； 如果A为多维数组，则沿第一个非单元素维进行处理，返回一个中间值数组。 5．求元素和函数sum 格式：s=sum(x) 功能：如果x为向量，则返回向量x的元素和；如果x为矩阵，则将矩阵每列当作向量处理，返回一个元素分别为各列和的行向量；如果A为多维数组，则沿第一个非单元素维进行计算，返回一个元素和数组。

6．求标准偏差函数std与方差函数var 对于向量 有两种标准差定义方法： 其中 为样本的元素个数。 格式一：s=std(x) 功能：x为向量，则返回用s1计算的标准偏差s。如果x是服从正态分布的随机样本，则s2为其方差的最佳无偏估计；如果x为矩阵，则返回矩阵每列标准差的行向量；如果x为多维数组，则沿x第一个非单元素维计算元素的标准偏差。 格式二：s=std(x, flag) 功能：如果flag=0，与s=std(x)一样 ；如果 flag =1，则返回用s2计算的标准差。

7．排序函数sort 格式一：A=sort(x) 功能：沿数组的不同维，以升序排列元素。元素可以为实数、复数和字符串。如果x是一个复数，其元素按其模的大小进行排列，如果模相等，则按其在区间[-pi, pi]上的相角进行排序。 格式二：[A ,index]=sort(x) 功能：同时返回一个下标数组index。 8．元素乘积函数prod 格式一：A=prod(x) 功能：如果x为向量，则计算其所有元素的乘积；如果x为矩阵，则每列作为向量处理，回一个每列元素积的行向量；如果x为多维数组，则沿第一个非单元素维进行处理，返回元素积数组。 格式二：A=prod(x,dim) 功能：沿dim指定维，返回元素积。

9．列元素累乘积函数cumprod( ) 格式一：A=cumprod(x) 功能：沿数组不同维，返回累乘积，返回值A与x大小一样，与元素全乘积不同，它只将x中相应元素与其之前的所有元素相乘。当x是向量，返回x的元素累计积向量；如果x为矩阵，返回一个与x大小相同的每列累乘积的矩阵；如果x为多维数组则沿第一个非单元素维计算累乘积。 格式二：A=cumprod(x, dim) 功 能：返回沿dim指足的维的元素的累乘积。 10．累计和cumsum函数 格式一：A=cumsum(x) 功 能：沿数组不同维，返回累乘和。当x是向量，返回x的元素累计和；如果x为短阵，返回一个与x大小相同按列累计和的矩阵；如果x为多维数组则沿第一个非单元素维计算累计和。 格式二：A=cumsum(x,dim) 功 能：沿dim指定的维计算元素的累计和。

3.3.2 协方差与相关系数 1．求协方差函数cov( ) 协方差函数定义为： 其中E表示数学期望； 格式一：C=cov(x) 功能：如果x为向量，则返回向量元素的方差；如果x为矩阵，每列产生一个方差向量，cov(x) 是一个协方差矩阵，diag(cov(x))为每列的方差向量，sqrt(diag(cov(x)))是一个标准差向量。 格式二：C=cov(x, y) 功能：返回x、y的协方差。x、y为长度相同的列向量。也可用C=cov([x,y])。

2．求相关系数函数corrcoef( ) 格式一：S=corrcoef(x) 功能：根据输入矩阵x，返回一个相关系数矩阵，相关系数S的矩阵与协方差矩阵 C=cov(x)有关，由下式确定： 格式二：S=corrcoef(x,y) 功能：返回列向量x和y的相关系数，也可用 S=corrcoef([x y])。

3.3.3 有限差分 1．求元素之差函数diff( ) 格式一：A=diff(x) 功能：计算x中相邻元素之间的差值或近似导数。如果x为向量，则返回一个比x少一个元素的向量，其元素值为[x(2)-x(1)，x(3)-x(2)，…，x(n)-x(n-1)]；如果x为矩阵，则返回一个列间差值的矩阵：[x(2:n,:) - x(1:n-1,:)]。 格式二：A=diff(x, n) 功能：使用diff函数递归 n 次，计算第n阶差值。例如，diff(x,2) = diff(diff(x))。 格式三：A=diff(x, n, dim) 功能：沿dim指定的维数计算第n阶差值。如果n大于或等于dim维的长度，则返回空数组。

2．求数值梯度函数gradient( ) 两变量函数F(x, y)的梯度定义为 对N个变量函数 F(x,y,z, … ) 其梯度为 梯度可看作指向F增加方向的向量集。 格式一：Fx=gradient(F) 功能：F为一向量，返回F的一维数值梯度，Fx与 一致，表示x方向的差分。 格式二：[Fx, Fy]=gradient(F) 功能：F为一矩阵，返回二维数值梯度的x和y分量。 Fx与 表示x(列)方向的差分 Fy与 一致，表示y(行)方向的差分。每个方向点间距离设为1。

3.4 函数分析与数值积分 3.4.1 函数在MATLAB中的表示与函数的绘图 1．函数的表示与计算 2．函数的绘制 (1) 单变量函数绘画命令fplot fplot( )函数的功能能确保在输出的图形中表示出所有的奇异点。 格式一：fplot(‘fun’, [xmin, xmax]) 功能：变量在[xmin, xmax]范围，绘制指定函数’fun’的图形。 格式二：fplot(‘fun’, [xmin, xmax], tol) 功能：在[xmin, xmax]范围，以给定的精度tol<1，绘制指定函数’fun’的图形，tol的缺省值为2e-3。

功能：当f=f(x)时，绘制函数f=f(x)在a < x < b范围函数曲线，并且所绘图上还自动进行 (2) 简易的函数绘图命令ezplot( ) 格式一：ezplot(f, [a,b]) 功能：当f=f(x)时，绘制函数f=f(x)在a < x < b范围函数曲线，并且所绘图上还自动进行 标注；当只输入函数文件名，而没有规定自变量的范围，其默认的自变量范围为 当f = f(x,y)时，则绘制a < x < b 、 a < y < b范围内f(x,y) = 0的函数曲线； 当只输入函数文件名，而没有规定自变量的范围，其默认的自变量范围为 格式二：ezplot(f, [xmin,xmax,ymin,ymax]) 功能：在xmin < x < xmax、ymin < y < ymax范围绘制f(x,y) = 0的曲线。 格式三：ezplot(x,y, [tmin,tmax]) 功能：在tmin < t < tmax范围绘制x = x(t)和y = y(t)的曲线。范围缺省是变量t的范围为 0 < t < 2*pi。 、

3.4.2函数的极点、零点分析 1．极值分析函数 (1) 单变量函数求极小值函数fminbnd( ) 格式：x=fminbnd (‘fun’,x1,x2) 功能：返回函数fun(x)在区间[x1, x2]内的局部极小值。 (2) 多变量函数求极小值函数fminsearch( ) fminsearch函数与fminbnd函数类似，但是它面向多变量函数。 格式：x=fminsearch(‘fun’,x0) 功能：返回x0附近，函数fun的局部极小化向量x。x0可以是标量、向量或矩阵。 2．单变量函数的零点分析 格式：x=fzero(‘funname’,x0) 功能：在x0附近，寻找函数funname的零点。funnamce为一个函数名的字符串，函数为单变量实值函数。Funnamce可以为函数句柄，也可以是inline对象。函数返回值的附近函数变号。如果x为两元素向量 , 则认为x0为区间，f(x0(1))的符号与f(x0(2))的符号相反，否则返回NaN。如果找到Inf、NaN，或复数值，则停止在查找区间内的搜索。

3.4.3函数的数值积分与微分 1．函数的数值积分 (1) 低阶数值积分函数quad 格式一：q=quad(‘fun’,a,b) 功能：采用自适应的Simpson积分方法，返回函数‘fun’ 在上限 a 和下限 b之间的数值积分。当给定一个输入值向量，‘fun’ 必须返回一个输出向量。函数‘fun’可以是函数名、函数句柄或字符串。 格式二：q=quad(‘fud’,a,b,tol) 功能：按指定绝对误差tol返回数值积分值, tol 缺省值为 1e-6。 (2) 高阶数值积分函数quadl

(3) 梯形面积法的积分函数trapz( ) 格式一：T=trapz(Y) 功能：以单位间隔，采用计算若干梯形面积的和来计算某函数的近似积分。如果Y为向量，计算Y的积分；如果Y是矩阵，得一个每列积分的行向量；如果Y为多维数组，则沿第一个非单元素维计算。 格式二：T=trapz(X,Y) 功能：用梯形积分法，依据X计算Y的积分。如果X为矢量，则Y必须是同大小的矢量；如果X是一列向量，并且数组Y第一非单元素维长度为length(X)，则在该维中计算。 (4) 双重积分函数dblquad MATLAB提供了一个求双重积分的函数dblquad，其基本调用格式为： 格式：Q=dblquad(fun,xmin,xmax,ymin,ymax,tol) 功能：按指定精度tol，对指定函数 f(x, y)在[xmin, xmax]范围和[ymin, ymax]范围进行双重积分。精度tol缺省时默认精度为1e-6。

(5) 不定积分的计算 对于函数不定积分的计算，可以采用定积分函数来求不定积分的数值解。方法是： 固定积分下限，用for循环，把积分上限逐步增加即可实现。 2．函数的数值微分 (1)通过计算数组中元素间的差分函数diff来粗略计算微分函数。 (2)对数据拟合后再利用polyder函数微分 为了实现对数值较准确的微分，最好先用最小二乘曲线拟合这种数据，然后对所 得到的多项式进行微分。

3.4.4 常微分方程的数值求解 任何高阶常微分方程都可以变换成一阶微分方程，即表示成右函数形式，这是利用龙格-库塔法求解微分方程的前提。 1．初始值的常微分方程求解 格式：[T,Y] = ode45(‘odefun’,tspan,y0) 功能：返回由文件‘odefun’所定义的具有初始条件为y0、时间t变化范围为[t0, tfinal] 的微分方程y‘ = f(t,y)的解，其中tspan = [t0,tfinal]。向量T中的每一列对应着矩阵Y的每一列。

2．边界条件的常微分方程求解 对于如下的微分方程： 用于边界条件的常微分方程求解问题: 函数bvp4c 格式：sol = bvp4c(‘odefun’,bcfun,solinit) 其中‘odefun’为常微分方程函数，bcfun为边界条件函数，solinit为求解的初始值。输出sol是一个结构，它有4个属性。 ·sil.x　bvp4v选择的网格节点 ·sil.y 网格点sil.x处y(x)的近似值 ·sil.yp　网格点sil.x处y’(x)的近似值 ·sil.parameters　未知参数的值。如果存在未知参数，则求解函数会自动求得。

3．偏微分方程的求解 考虑如下的偏微分方程： 其中有t0<t<tf ，a<x<b,，m=0、1或2。如果m>0，则必须a>0，求满足初始条件 u(x, t0)=u0(x)边界条件如下的解： MATLAB提供了dpepe函数来求解该问题的数值解。其基本调用格式为： sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)