1、第五讲第五讲MATLABMATLAB数值计算数值计算5.1 5.1 特殊矩阵特殊矩阵5.2 5.2 矩阵分析矩阵分析5.3 5.3 矩阵分解与线性方程组求解矩阵分解与线性方程组求解5.4 5.4 数据处理与多项式计算数据处理与多项式计算5.5 5.5 傅立叶分析傅立叶分析5.6 5.6 数值微积分数值微积分5.7 5.7 常微分方程的数值求解常微分方程的数值求解5.8 5.8 非线性方程的数值求解非线性方程的数值求解5.9 5.9 稀疏矩阵稀疏矩阵5.1 特殊矩阵5.1.1 5.1.1 对角阵与三角阵对角阵与三角阵1.1.矩阵的对角元素矩阵的对角元素(1)提取矩阵的对角线元素提取矩阵的对角线元
2、素 设设A为为mn矩阵，矩阵，diag(A)函数用于提取矩阵函数用于提取矩阵A主对主对角线元素产生一个具有角线元素产生一个具有min(m,n)个元素的列向量。个元素的列向量。(2)构造对角矩阵构造对角矩阵 设设V为具有为具有m个元素的向量，个元素的向量，diag(V)将产生一个将产生一个mm对角矩阵，其主对角线元素即为向量对角矩阵，其主对角线元素即为向量V的元素。的元素。例例5.1 先建立先建立55矩阵矩阵A，然后将，然后将A的第的第1行元素乘行元素乘以以1，第，第2行乘以行乘以2，第，第5行乘以行乘以5。命令如下：A=17,0,1,0,15;.23,5,7,14,16;.4,0,13,0,2
3、2;.10,12,19,21,3;.11,18,25,2,19D=diag(1,2,3,4,5)D*A 2.2.矩阵的三角阵矩阵的三角阵 (1)下三角矩阵下三角矩阵 求矩阵求矩阵A的下三角阵的下三角阵:tril(A)。(2)上三角矩阵上三角矩阵 提取矩阵提取矩阵A的上三角矩阵的上三角矩阵:triu(A)5.1.2 5.1.2 特殊矩阵的生成特殊矩阵的生成 1.1.魔方矩阵魔方矩阵 函数函数magic(n):生成一个生成一个n阶魔方阵。阶魔方阵。magic(3)ans=8 1 6 3 5 7 4 9 2例例5.2 将将101125等等25个数个数填入一个填入一个5行行5列的表格中，列的表格中，使
4、其每行每列及对角线的使其每行每列及对角线的和均为和均为565。B=100+magic(5)A=vander(1:5)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 12.范得蒙矩阵范得蒙矩阵(Vandermonde)函数函数 A=vander(V)生成以向量生成以向量V为基础向为基础向量的范得蒙矩阵。量的范得蒙矩阵。)()()(jniViA 3.希尔伯特矩阵希尔伯特矩阵(Hilbert)生成希尔伯特矩阵的函数生成希尔伯特矩阵的函数:H=hilb(n)hilb(3)ans=1.0000 0.5000 0.3333 0.5
5、000 0.3333 0.2500 0.3333 0.2500 0.2000求求n阶希尔伯特矩阵的逆的函数阶希尔伯特矩阵的逆的函数:invhilb(n)11),(jijiH 4.托普利兹矩阵托普利兹矩阵(Toeplitz)生成托普利兹矩阵的函数生成托普利兹矩阵的函数:T=toeplitz(c,r)，它生成一个以，它生成一个以c为第为第1列，列，r为第为第1行的托普利兹矩阵。这里行的托普利兹矩阵。这里c,r均为向量，二均为向量，二者不必等长。者不必等长。c=1 2 3 4 5;r=1.5 2.5 3.5 4.5 5.5;toeplitz(c,r)Warning:First element of
6、input column does not match first element of input row.Column wins diagonal conflict.ans=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.伴随矩阵伴随矩阵 生成伴随矩阵的函数是：生成伴随矩阵的函数是：A=co
7、mpan(P)，生成多项，生成多项式式P的伴随矩阵。的伴随矩阵。P是一个多项式的系数向量，高次是一个多项式的系数向量，高次幂系数排在前，低次幂排在后。伴随矩阵的特征值幂系数排在前，低次幂排在后。伴随矩阵的特征值是多项式的根。是多项式的根。u=1 0 -7 6A=compan(u)A=0 7 -6 1 0 0 0 1 067)3)(2)(1(3xxxxxeig(compan(u)ans=-3.0000 2.0000 1.00006.帕斯卡矩阵帕斯卡矩阵(Pascal)函数函数A=pascal(n)生成一个生成一个n阶的帕斯卡矩阵阶的帕斯卡矩阵。pascal(5)ans=1 1 1 1 1 1 2
8、 3 4 5 1 3 6 10 15 1 4 10 20 35 1 5 15 35 70),1()1,(),(,1)1,(,1),1(jiajiajiAiAjA5.2 5.2 矩阵分析矩阵分析5.2.1 5.2.1 矩阵结构变换矩阵结构变换1.矩阵的转置矩阵的转置转置运算符是单撇号转置运算符是单撇号()。2.矩阵的旋转矩阵的旋转矩阵的旋转利用函数矩阵的旋转利用函数rot90(A,k)，功能是将矩阵，功能是将矩阵A旋旋转转90的的k倍，当倍，当k为为1时可省略。时可省略。3.矩阵的左右翻转矩阵的左右翻转对矩阵对矩阵A实施左右翻转的函数是实施左右翻转的函数是fliplr(A)。4.矩阵的上下翻转矩
9、阵的上下翻转对矩阵对矩阵A实施上下翻转的函数是实施上下翻转的函数是flipud(A)。5.2.2 矩阵的逆与伪逆矩阵的逆与伪逆1.矩阵的逆矩阵的逆 求方阵求方阵A的逆可调用函数的逆可调用函数inv(A)。例例5.4 用求逆矩阵的方法解线性方程组。用求逆矩阵的方法解线性方程组。A=1,2,3;1,4,9;1,8,27;b=5,2,6;x=inv(A)*b一般情况下，用左除比求矩阵的逆的方法更一般情况下，用左除比求矩阵的逆的方法更有效，即有效，即x=Ab。2.矩阵的伪逆矩阵的伪逆对奇异方阵和长方阵，求矩阵伪逆的函数是对奇异方阵和长方阵，求矩阵伪逆的函数是pinv(A)。例例5.5 求求A的伪逆，并
10、将结果送的伪逆，并将结果送B。A=3,1,1,1;1,3,1,1;1,1,3,1;B=pinv(A)例例5.6 求矩阵求矩阵A的伪逆。的伪逆。A=0,0,0;0,1,0;0,0,1;pinv(A)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);%用方程组的右端向量置换用方程组的右
11、端向量置换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 5.2.4 矩阵的秩矩阵的秩求矩阵秩的函数求矩阵秩的函数:rank(A)。例如，求例例如，求例5.7中方程组系数
12、矩阵中方程组系数矩阵D的秩，命令是：的秩，命令是：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;r=rank(D)r=4说明说明D是一个满秩矩阵。是一个满秩矩阵。5.2.5 向量和矩阵的范数向量和矩阵的范数1.计算向量计算向量3种常用范数的函数种常用范数的函数(1)norm(V)或或norm(V,2)计算向量计算向量V的的2范数范数 (sum(abs(V).2)(1/2)(2)norm(V,1)计算向量计算向量V的的1范数范数(sum(abs(V)(3)norm(V,inf)计算向量计算向量V的的范数范数(max(abs(V)例例5.8 已知已知V，求，求V的的3
13、种范数。种范数。V=-1,1/2,1;v1=norm(V,1)%求求V的的1范数范数v2=norm(V)%求求V的的2范数范数vinf=norm(V,inf)%求求范数范数 2.矩阵的范数及其计算函数矩阵的范数及其计算函数MATLAB中提供了求中提供了求3种矩阵范数的函数，其函数调用格式与求种矩阵范数的函数，其函数调用格式与求向量的范数的函数完全相同。向量的范数的函数完全相同。例例5.9 求矩阵求矩阵A的三种范数。的三种范数。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的的
14、1-范数范数(列和范数列和范数):max(sum(abs(A)a2=norm(A)%求求A的的2-范数范数(谱范数谱范数):max(svd(A)ainf=norm(A,inf)%求求A的的-范数范数(行和范数行和范数):max(sum(abs(A)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的的 范数下的条件数
15、范数下的条件数例例5.10 求矩阵求矩阵X的三种条件数。的三种条件数。A=2,2,3;4,5,-6;7,8,9;C1=cond(A,1)C2=cond(A)C3=cond(A,inf)2.矩阵的迹矩阵的迹求矩阵的迹的函数求矩阵的迹的函数:trace(A)。例如，例如，X=2 2 3;4 5-6;7 8 9;trace(X)ans=165.2.7 矩阵的特征值与特征向量矩阵的特征值与特征向量计算矩阵计算矩阵A的特征值和特征向量的函数是的特征值和特征向量的函数是eig(A)，常用的调用格式有常用的调用格式有3种：种：(1)E=eig(A)求矩阵求矩阵A的全部特征值，构成向量的全部特征值，构成向量E
16、。(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,1E=eig(A)V,
17、D=eig(A)V,D=eig(A,nobalance)例例5.12 用求特征值的方法解方程。用求特征值的方法解方程。p=3,-7,0,5,2,-18;%A=compan(p);%A的伴随矩阵的伴随矩阵x1=eig(A)%求求A的特征值的特征值x2=roots(p)%直接多项式直接多项式p的零点的零点两种方法求得的方程的根是完全一致的，实际上，两种方法求得的方程的根是完全一致的，实际上，roots函数正是应用求伴随矩阵的特征值的方法函数正是应用求伴随矩阵的特征值的方法来求方程的根。来求方程的根。5.2.8 MATLAB在三维向量中的应用在三维向量中的应用1.向量共线或共面的判断向量共线或共面的
18、判断例例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)2.向量方向余弦的计算向量方向余弦的计算例例5.14 设向量设向量V=(5，-3，2)，求，求V的方向余弦。的方向余弦。建立一个函数文件建立一个函数文件direct.m：function f=direct(v)r=norm(v);if r=0 f=0else f=v(1)/r
19、,v(2)/r,v(3)/r;end在在MATLAB命令窗口，输入命令：命令窗口，输入命令：V=5,-3,2;F=direct(v)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.向量的向量积向量
20、的向量积例例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=65.3 5.3 矩阵分解矩阵分解1.实对称矩阵的实
21、对称矩阵的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 求下列二次型的标准形式及变换矩阵。求下列二次型的标准形式及变换矩阵。命令如下
22、：命令如下：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。3.矩阵的矩阵的QR分解分解对矩阵对矩阵A进行进行QR分解的函数是分解的函数是Q,R=qr(A)，根据方
23、阵根据方阵A，求一个，求一个正交矩阵正交矩阵Q和一个和一个上三上三角矩阵角矩阵R，使，使A=Q*R。例如，对矩阵例如，对矩阵A进行进行QR分解的命令是：分解的命令是：A=2,1,-2;1,2,1;2,5,3;Q,R=qr(A)5.4 5.4 数据处理与多项式计算数据处理与多项式计算5.4.1 5.4.1 数据统计与分析数据统计与分析1.求矩阵最大和最小元素求矩阵最大和最小元素(1)求向量的最大最小元素求向量的最大最小元素y=max(X)返回向量返回向量X的最大元素存入的最大元素存入y。y,I=max(X)返回向量返回向量X的最大元素存入的最大元素存入y，最大元素的，最大元素的序号存入序号存入I
24、。(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行上的最大元素。行上的最大元素。(
25、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的每行及每列的最大和最小元素，的每行及每列的最大和最小元素，并求整个矩
26、阵的最大和最小元。并求整个矩阵的最大和最小元。命令如下：命令如下：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)%求整个矩阵的最小元素求整个矩阵的最小元素2.求矩阵的平均值和中值求矩阵的平均值和中值求矩阵和向量元素的平均值的函数是求矩阵和向量元素的平均值的函数是mean，求，求中值的函数是中值的函数是med
27、ian。它们的调用方法和。它们的调用方法和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函
28、数能方便函数能方便地求得向量和矩阵元素的累加和与累乘积向量，地求得向量和矩阵元素的累加和与累乘积向量，函数的用法和函数的用法和sum及及prod相同相同例例5.27求向量求向量X=(1！，！，2！，！，3！，！，10！)。命令如下：命令如下：X=cumprod(1:10)5.标准方差标准方差MATLAB中，提供了计算数据序列的标准方差的函中，提供了计算数据序列的标准方差的函数数std。对于向量对于向量X，std(X)返回一个标准方差。对返回一个标准方差。对于矩阵于矩阵A，std(A)返回一个行向量，它的各个元素返回一个行向量，它的各个元素便是矩阵便是矩阵A各列或各行的标准方差。各列或各行的标准
29、方差。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,mo
30、de)其中其中dim指明对指明对A的列还是行进行排序，若的列还是行进行排序，若dim=1，则 按 列 排 序，若则 按 列 排 序，若 d i m=2，则 按 行 排 序。，则 按 行 排 序。mode=ascend按升序排序按升序排序(缺省缺省),mode=descend按降序排序，按降序排序，Y是排序后的矩阵，而是排序后的矩阵，而I记录记录Y中的中的元素在元素在A中位置。中位置。例例5.28 对矩阵做各种排序。对矩阵做各种排序。命令如下：命令如下：A=1,-8,5;4,12,6;13,7,-13;sort(A)%对对A的每列按升序排序的每列按升序排序sort(A,2,descend)%对对
31、A的每行按降序排序的每行按降序排序X,I=sort(A)%对对A按列排序，并将每个元素按列排序，并将每个元素所在行号送矩阵所在行号送矩阵I5.4.2 数值插值数值插值1.一维数值插值一维数值插值interp1函数调用格式为：函数调用格式为：Y1=interp1(X,Y,X1,method)函数根据函数根据X、Y的值，计算函数在的值，计算函数在X1处的值。处的值。X、Y是两个等长的已知向量，分别描述采样点和样是两个等长的已知向量，分别描述采样点和样本值，本值，X1是一个向量或标量，描述欲插值的点，是一个向量或标量，描述欲插值的点，Y1是一个与是一个与X1等长的插值结果。等长的插值结果。metho
32、d是插值是插值方法，允许的取值有方法，允许的取值有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)i
33、nterp1(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)M A T L A B 中 有 一 个 专 门 的 三 次 样 条 插 值 函 数中 有 一 个 专 门 的 三 次 样 条 插 值 函 数Y1=spline(X,Y,X1)，其功能及使用方法与函数，其功能及使用方
34、法与函数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.2015,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=inter
35、p1(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是与参数采
36、样点对应的采样变量的样本是与参数采样点对应的采样变量的样本值，值，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)%对函数在对函数在
37、(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)点进行插值点进行插值3.三维数值插值三维数值插值对三维函数插值的函数是对三维函数插值的函数是interp3，其使用方法和，其使用方法和interp2相同。其调用格式为：相同。其调用格式为：W1=interp3(X,Y,Z,W,X1,Y1
38、,Z1,method)函数返回三维插值结果。其中函数返回三维插值结果。其中X、Y、Z是三个向是三个向量，分别描述三个参数的采样点，量，分别描述三个参数的采样点，W是与参数是与参数采样点对应的采样变量的样本值，采样点对应的采样变量的样本值，X1、Y1、Z1是三个向量或标量，描述欲插值的点。是三个向量或标量，描述欲插值的点。method是插值方法，可选，其缺省值是是插值方法，可选，其缺省值是 line。method的取值与一、二维插值函数相的取值与一、二维插值函数相同。同。5.4.3 曲线拟合曲线拟合(一一)多项式曲线拟合多项式曲线拟合MATLAB中，提供了解决使用最小二乘法进行曲中，提供了解决使
39、用最小二乘法进行曲线拟合的函数。调用格式为：线拟合的函数。调用格式为：P,S=polyfit(X,Y,m)函数根据采样点函数根据采样点X和采样点函数值和采样点函数值Y，产生一个，产生一个m次多项式次多项式P及其在采样点的误差向量及其在采样点的误差向量S。其中其中X、Y是两个等长的向量，是两个等长的向量，P是一个长度为是一个长度为m+1的向量。的向量。P(1)*xm+P(2)*x(m-1)+P(m)*x+p(m+1)例例5.32 用一个用一个5次多项式在区间次多项式在区间0，2内逼近函数内逼近函数sin(x)。命令如下：命令如下：X=linspace(0,2*pi,50);Y=sin(X);P,
40、S=polyfit(X,Y,5)%得到得到5次多项式的系数和误差次多项式的系数和误差plot(X,Y,k*,X,polyval(P,X),k-)（二）使用指定函数进行曲线拟合（二）使用指定函数进行曲线拟合x=lsqnonlin(fun,x0)x返回拟合参数，返回拟合参数，fun为拟合曲线对应的函数，为拟合曲线对应的函数，x0为初始值为初始值例：根据所给试验数据，用双曲线模型例：根据所给试验数据，用双曲线模型拟合，确定模型参数拟合，确定模型参数0 xyyabx0,0.05,0.1,0.2,0.4,0.6,0.8,1.21.335,1.253,1.180,1.058,0.887,0.803,0.7
41、52,0.685xyfunction=curvefit()data=0 1.335;0.05 1.253;0.1 1.180;0.2 1.058;0.4 0.887;0.6 0.803;0.8 0.752;1.2 0.685;x=data(:,1);y=data(:,2);plot(x,y,g-);hold on;x0=1 0;p=lsqnonlin(curve_model,x0);y1=y(1);y=y1-x./(p(1)+p(2)*x);plot(x,y,bo);%定义curve_model函数function f=curve_model(p,data)data=0 1.335;0.05
42、1.253;0.1 1.180;0.2 1.058;0.4 0.887;0.6 0.803;0.8 0.752;1.2 0.685;x=data(:,1);y=data(:,2);y1=y(1);z=y1-x./(p(1)+p(2)*x);f=z-y;5.4.4 多项式计算多项式计算1.多项式的建立多项式的建立已知一个多项式的全部根已知一个多项式的全部根X求多项式系数的函数是求多项式系数的函数是poly(X)，该函数返回以，该函数返回以X为全部根的一个多项式为全部根的一个多项式P，当当X是一个长度为是一个长度为m的向量时，的向量时，P是一个长度为是一个长度为m+1的向量。的向量。2.多项式求根
43、多项式求根求多项式求多项式p(x)的根的函数是的根的函数是roots(P)，这里，这里，P是是p(x)的系数向量，该函数返回方程的系数向量，该函数返回方程p(x)=0的全部根的全部根(含含重根，复根重根，复根)。3.多项式求值多项式求值求多项式求多项式p(x)在某点或某些点的函数值的函数是在某点或某些点的函数值的函数是polyval(P,x)。若。若x为一数值，则求多项式在该点的为一数值，则求多项式在该点的值；若值；若x为向量或矩阵，则对向量或矩阵中的每个为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。元素求其多项式的值。例例5.33 已知一个多项式，已知一个多项式，(1)计算计算f(
44、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,5X=roots(P)%求方程求方程f(x)=0的根的根G=poly(X)%求多项式求多项式g(x)X0=5,7.8,9.6,12.3f=polyval(P,X0)%求多项式求多项式f(x)在给定点的值在给定点的值多项式求值还有一个函数是多项式求值还有一个函数是polyvalm，其调用格式，其调用格式与与polyval相同，
45、但含义不同。相同，但含义不同。polyvalm函数要求函数要求x为方阵，它以方为方阵，它以方阵为自变量求多项式的值。阵为自变量求多项式的值。4.多项式的四则运算多项式的四则运算(1)多项式的加减法多项式的加减法(2)多项式的乘法多项式的乘法函数函数conv(P1,P2)用于求多项式用于求多项式P1和和P2的乘积。的乘积。(3)多项式的除法多项式的除法函数函数Q,r=deconv(P1,P2)用于对多项式用于对多项式P1和和P2作除作除法运算。其中法运算。其中Q返回多项式返回多项式P1除以除以P2的商式，的商式，r返回返回P1除以除以P2的余式。这里，的余式。这里，Q和和r仍是多项式仍是多项式系
46、数向量。系数向量。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.多项式的导函数
47、多项式的导函数对多项式求导数的函数是：对多项式求导数的函数是：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中用于求最小值的函数是：中用于求最小值的函数是：fminbnd 求一元
48、函数在一确定区间内的极值。求一元函数在一确定区间内的极值。fminsearch 求多元函数的极值。求多元函数的极值。vMATLAB没有专门提供求函数最大值点的函数，没有专门提供求函数最大值点的函数，但只要注意到但只要注意到-f(x)在区间在区间(a,b)上的最小值点就是上的最小值点就是f(x)在在(a,b)的最大值点，所以的最大值点，所以fmin(-f,a,b)返回函返回函数数f(x)在区间在区间(a,b)上的最大值。上的最大值。x,feval,exitflag,output=fminbnd(fun,x1,x2,options)x为函数为函数fun在区间在区间x1x0，函数收敛到函数收敛到x；
49、exitflag=0，函数计算次数已达最大迭代次数；，函数计算次数已达最大迭代次数；exitflag0，函数在计算区间内不收敛。函数在计算区间内不收敛。output为包含最优化信息的结构，为包含最优化信息的结构，output.algorithm，output.funccount,output.iterations,output.message分别为分别为所用算法、函数计算次数，迭代次数，退出信息所用算法、函数计算次数，迭代次数，退出信息options为可选项，如：为可选项，如：Display为显示层次：为显示层次：off时不显示时不显示输出内容；输出内容；iter时显示迭代过程，时显示迭代过程
50、，final时显示输出内容；时显示输出内容；notify时函数不收敛则显示输出（默认选择）。时函数不收敛则显示输出（默认选择）。MaxFunEvals为函数值最大误差。为函数值最大误差。MaxIter为迭代最大步为迭代最大步数数(默认默认500)。TolX为为x误差范围（默认误差范围（默认1.0e-4)例例5.36 求函数求函数 在区间在区间(0，2)上的极值。上的极值。3()25f xxxf=inline(x.3-2*x-5);%通过内联函数建立函数通过内联函数建立函数f x1=fminbnd(f,0,2)%求函数求函数f在区间在区间(0,2)上的极值上的极值x2=fminbnd(f,0,2
