第4章__MATLAB_数值计算(1)20140822

1、123u数据处理与多项式计算数据处理与多项式计算u数值微积分数值微积分u离散傅立叶变换离散傅立叶变换u线性方程组求解线性方程组求解u非线性方程与最优化问题求解非线性方程与最优化问题求解u常微分方程的数值求解常微分方程的数值求解u稀疏矩阵稀疏矩阵4u数据统计与分析数据统计与分析 矩阵平均值与均值矩阵平均值与均值 标准差标准差 相关系数相关系数 矩阵最大和最小元素矩阵最大和最小元素 元素求和与求积元素求和与求积 元素排序元素排序u数据插值数据插值u曲线拟合曲线拟合u多项式计算多项式计算5u平均值平均值 算术平均值算术平均值u中值中值 数据序列中数据序列中其值的大小恰好处在中间的元素其值的大小恰好处

2、在中间的元素 例：奇数个数据序列例：奇数个数据序列-2，5，7，9，12的中值的中值 偶数个数据序列偶数个数据序列-2，5，6，7，9，12的中值的中值u求矩阵和向量元素的平均值和中值求矩阵和向量元素的平均值和中值 平均值的函数是平均值的函数是mean 中值的函数是中值的函数是median6u设设X是一个向量，是一个向量，A是一个矩阵，是一个矩阵，mean和和median两个函数的调用格式两个函数的调用格式为为 mean(X)：返回向量：返回向量X的算术平均值的算术平均值 median(X)：返回向量：返回向量X的中值。的中值。 mean(A)：返回一个行向量，其第：返回一个行向量，其第i个元

3、素是个元素是A的第的第i列的算术平均值。列的算术平均值。 median(A)：返回一个行向量，其第：返回一个行向量，其第i个元素是个元素是A的第的第i列的中值。列的中值。 mean(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于mean(A)；当；当dim为为2时，返回时，返回一个列向量，其第一个列向量，其第i个元素是个元素是A的第的第i行的算术平均值。行的算术平均值。 median(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于median(A)；当；当dim为为2时，返时，返回一个列向量，其第回一个列向量，其第i个元素是个元素是A的第的第i行的中值。行的

4、中值。u例：例：A=4,9,6;1,8,3; 7,2,9 mean(A) median(A) mean(A,2) median(A,2)7u数据序列的标准方差的计算公式数据序列的标准方差的计算公式u在在MATLAB中，提供了计算数据序列的标准方差的函数中，提供了计算数据序列的标准方差的函数std。 对于向量对于向量X，std(X)返回一个标准方差。返回一个标准方差。 对于矩阵对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵返回一个行向量，它的各个元素便是矩阵A各列或各列或各行的标准方差。各行的标准方差。21122111()11()1NiiNiiNiiSxxNSxxNxxN8ustd

5、函数的一般调用格式为：函数的一般调用格式为：Y=std(A,flag,dim)uDim的取值的取值=1时，求各列元素的标准方差时，求各列元素的标准方差=2时，求各行元素的标准方差时，求各行元素的标准方差uFlag的取值的取值=0时，按时，按S1所列公式计算标准方差所列公式计算标准方差=1时，按时，按S2所列公式计算标准方差所列公式计算标准方差缺省缺省flag=0，dim=1u例例6.4 对二维矩阵对二维矩阵x，从不同维方向求出其标准方差。，从不同维方向求出其标准方差。x=4,5,6;1,4,8 % x=4,5,6;1,4,8 % 产生一个二维矩阵产生一个二维矩阵x xy1=std(x,0,1)

6、y1=std(x,0,1)y2=std(x,1,1)y2=std(x,1,1)y3=std(x,0,2)y3=std(x,0,2)y4=std(x,1,2)y4=std(x,1,2)9u相关系数的计算公式相关系数的计算公式uMATLAB提供了提供了corrcoef函数，可以求出数据的相关系数矩函数，可以求出数据的相关系数矩阵。阵。ucorrcoef函数的调用格式为：函数的调用格式为： corrcoef(X)：返回从矩阵：返回从矩阵X形成的一个相关系数矩阵。此相关系数形成的一个相关系数矩阵。此相关系数矩阵的大小与矩阵矩阵的大小与矩阵X一样。它把矩阵一样。它把矩阵X的每列作为一个变量，然后求的每列

7、作为一个变量，然后求它们的相关系数它们的相关系数 corrcoef(X,Y)：在这里，：在这里，X,Y是向量，它们与是向量，它们与corrcoef(X,Y)的作用的作用一样。一样。22()()()()xixyiyrxixyiy10u例例 生成满足正态分布的生成满足正态分布的100005随机矩阵，然后求各列随机矩阵，然后求各列元素的均值和标准方差，再求这元素的均值和标准方差，再求这5列随机数据的相关系数列随机数据的相关系数矩阵。矩阵。 X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X)11uMATLAB提供的求数据序列的最大值和最小值的函数分别为

8、提供的求数据序列的最大值和最小值的函数分别为max和和min，两个函数的调用格式和操作过程类似。，两个函数的调用格式和操作过程类似。u求向量和矩阵的最大值求向量和矩阵的最大值 y=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，如果，如果X中包含复数元素，则按模中包含复数元素，则按模取最大值。取最大值。 y,I=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，最大值的序号存入，最大值的序号存入I，如果，如果X中中包含复数元素，则按模取最大值。包含复数元素，则按模取最大值。 max(A)：返回一个行向量，向量的第：返回一个行向量，向量的第i个元素是矩阵个元素是矩阵A的

9、第的第i列上的最大值列上的最大值 Y,U=max(A)：返回行向量：返回行向量Y和和U，Y向量记录向量记录A的每列的最大值，的每列的最大值，U向量向量记录每列最大值的行号。记录每列最大值的行号。 max(A,dim)：dim取取1或或2。dim取取1时，该函数和时，该函数和max(A)完全相同；完全相同；dim取取2时，该函数返回一个列向量，其第时，该函数返回一个列向量，其第i个元素是个元素是A矩阵的第矩阵的第i行上的最大值行上的最大值u求向量和矩阵的最小值求向量和矩阵的最小值 求向量求向量X的最小值的函数是的最小值的函数是min(X)，用法和，用法和max(X)完全相同。完全相同。12u例例

10、 分别求矩阵分别求矩阵A中各列和各行元素中的最大值，中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。并求整个矩阵的最大值和最小值。 A=13,-56,78;25,63,-235;78,25,563;1,0,-1;A=13,-56,78;25,63,-235;78,25,563;1,0,-1; max(A,2) %max(A,2) %求每行最大元素求每行最大元素 min(A,2) %min(A,2) %求每行最小元素求每行最小元素 max(A) %max(A) %求每列最大元素求每列最大元素 min(A) %min(A) %求每列最小元素求每列最小元素 max(max(A) %max(

11、max(A) %求整个矩阵的最大元素。也可使用命令：求整个矩阵的最大元素。也可使用命令：max(A(:)max(A(:) min(min(A) %min(min(A) %求整个矩阵的最小元素。也可使用命令：求整个矩阵的最小元素。也可使用命令：min(A(:)min(A(:)1 35 67 82 56 32 3 57 82 55 6 3101A13u两个向量或矩阵对应元素的比较两个向量或矩阵对应元素的比较 U=max(A,B)：A,B是两个同型的向量或矩阵，结果是两个同型的向量或矩阵，结果U是与是与A,B同同型的向量或矩阵，型的向量或矩阵，U的每个元素等于的每个元素等于A,B对应元素的较大者。对

12、应元素的较大者。 U=max(A,n)：n是一个标量，结果是一个标量，结果U是与是与A同型的向量或矩阵，同型的向量或矩阵，U的每个元素等于的每个元素等于A对应元素和对应元素和n中的较大者。中的较大者。u例例 求两个求两个23矩阵矩阵X, Y所有同一位置上的较大元素构成的新矩阵所有同一位置上的较大元素构成的新矩阵P X=4,5,6;1,4,8 Y=1,7,5;4,5,7 P=max(X, Y )14u数据序列求和与求积的函数是数据序列求和与求积的函数是sum和和prod，其使用方法类似。，其使用方法类似。u设设X是一个向量，是一个向量，A是一个矩阵，函数的调用格式为：是一个矩阵，函数的调用格式为

13、：sum(X)：返回向量：返回向量X各元素的和各元素的和prod(X)：返回向量：返回向量X各元素的乘积。各元素的乘积。sum(A)：返回一个行向量，其第：返回一个行向量，其第i个元素是个元素是A的第的第i列的元素和。列的元素和。prod(A)：返回一个行向量，其第：返回一个行向量，其第i个元素是个元素是A的第的第i列的元素乘积。列的元素乘积。sum(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于sum(A)；当；当dim为为2时，返回一个列向量时，返回一个列向量，其第，其第i个元素是个元素是A的第的第i行的各元素之和。行的各元素之和。prod(A,dim)：当：当dim为为

14、1时，该函数等同于时，该函数等同于prod(A)；当；当dim为为2时，返回一个列向时，返回一个列向量，其第量，其第i个元素是个元素是A的第的第i行的各元素乘积。行的各元素乘积。例例 求矩阵求矩阵A的每行元素的乘积和全部元素的乘积。的每行元素的乘积和全部元素的乘积。A=1,2,3,4;5,6,7,8;9,10,11,12;A=1,2,3,4;5,6,7,8;9,10,11,12;S=prod(A,2)S=prod(A,2)prod(S) %prod(S) %求求A A的全部元素的乘积。也可以使用命令的全部元素的乘积。也可以使用命令prod(A(:)prod(A(:)15u在在MATLAB中，使

15、用中，使用cumsum和和cumprod函数能方便地求得向量和矩函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的调用格式为：阵元素的累加和与累乘积向量，函数的调用格式为： cumsum(X)：返回向量：返回向量X累加和向量。累加和向量。 cumprod(X)：返回向量：返回向量X累乘积向量。累乘积向量。 cumsum(A)：返回一个矩阵，其第：返回一个矩阵，其第i列是列是A的第的第i列的累加和向量。列的累加和向量。 cumprod(A)：返回一个矩阵，其第：返回一个矩阵，其第i列是列是A的第的第i列的累乘积向量。列的累乘积向量。 cumsum(A,dim)：当：当dim为为1时，该函

16、数等同于时，该函数等同于cumsum(A)；当；当dim为为2时时，返回一个矩阵，其第，返回一个矩阵，其第i行是行是A的第的第i行的累加和向量。行的累加和向量。 cumprod(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于cumprod(A)；当；当dim为为2时，返回一个向量，其第时，返回一个向量，其第i行是行是A的第的第i行的累乘积向量。行的累乘积向量。u例：求向量例：求向量X=(1X=(1！，！，2 2！，！，3 3！，！，1010！) )。 X=cumprod(1:10) X=cumprod(1:10) 16uMATLAB中对向量中对向量X是排序函数是是排序函数是s

17、ort(X)，函数返回一个对，函数返回一个对X中的元素按升序中的元素按升序排列的新向量。排列的新向量。usort函数也可以对矩阵函数也可以对矩阵A的各列或各行重新排序，其调用格式为：的各列或各行重新排序，其调用格式为： Y,I=sort(A,dim)dim指明对指明对A的列还是行进行排序。若的列还是行进行排序。若dim=1，则按列排；若，则按列排；若dim=2，则按行排，则按行排Y是排序后的矩阵是排序后的矩阵I记录记录Y中的元素在中的元素在A中位置。中位置。u例例 对矩阵对矩阵A A做各种排序做各种排序A=1,-8,5;4,12,6;13,7,-13;A=1,-8,5;4,12,6;13,7,

18、-13;sort(A) %sort(A) %对对A A的每列按升序排序的每列按升序排序-sort(-A,2) %-sort(-A,2) %对对A A的每行按降序排序的每行按降序排序X,I=sort(A) %X,I=sort(A) %对对A A按列排序，并将每个元素所在行号送矩阵按列排序，并将每个元素所在行号送矩阵I I137136124581A17u为什么要做数据插值为什么要做数据插值? ? 自变量自变量x x与因变量与因变量y=f(x)y=f(x)的关系式有时不能直接写出表达式的关系式有时不能直接写出表达式，而，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点只能得到函数在若干个点

19、的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。之外的函数值时，需要估计函数值在该点的值。 由由有限个已知数据点有限个已知数据点，构造一个解析表达式构造一个解析表达式，由此计算数据点之，由此计算数据点之间的函数值，称之为插值。间的函数值，称之为插值。 01234567-1-0.8-0.6-0.4-0.200.20.40.60.8118u数据插值的分类数据插值的分类 根据插值函数的自变量个数根据插值函数的自变量个数一维一维二维二维多维多维 根据插值函数的选取根据插值函数的选取线性插值线性插值临近点插值临近点插值多项式插值多项式插值样条插值样条插值19u一维数据插值一维

20、数据插值是对一元数据点是对一元数据点(xi,yi)(xi,yi)进行插值进行插值 u在在MATLAB中，实现一维数据插值的函数是中，实现一维数据插值的函数是interp1，其调用格式为：，其调用格式为： Y1=interp1(X,Y,X1,method)X,Y是两个等长的已知向量，分别描述采样点和样本值，是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，是一个向量或标量，描述欲插值的点，Y1是一个与是一个与X1等长的插值结果。等长的插值结果。method是插值方法，允许的取值有是插值方法，允许的取值有linear、nearest、cubic、spline。u选

21、择插值方法时主要考虑因素：选择插值方法时主要考虑因素： 运算时间、占用计算机内存和插值的光滑程运算时间、占用计算机内存和插值的光滑程度。度。 运算时间运算时间 占用计算机内存占用计算机内存 光滑程度光滑程度 临近点插值临近点插值nearest ： 快快 少少 差差线性插值线性插值linear ： 稍长稍长 较多较多 稍好稍好三次样条插值三次样条插值splin： 最长最长 较多较多 最好最好立方插值立方插值cubic ： 较长较长 多多 较好较好20u例例: :已知数据：已知数据：u求当求当xi=0.25xi=0.25时的时的yiyi的值。的值。 x=0:.1:1;x=0:.1:1; y=.3

22、.5 1 1.4 1.6 1 .6 .4 .8 1.5 2;y=.3 .5 1 1.4 1.6 1 .6 .4 .8 1.5 2; yi0=interp1(x,y,0.25,linear)yi0=interp1(x,y,0.25,linear)u比较几种插值方法比较几种插值方法 xi=0:.02:1;xi=0:.02:1; yi=interp1(x,y,xi,linear);yi=interp1(x,y,xi,linear); zi=interp1(x,y,xi,spline);zi=interp1(x,y,xi,spline); wi=interp1(x,y,xi,cubic);wi=int

23、erp1(x,y,xi,cubic); plot(x,y,o,xi,yi,r+, xi,zi,gplot(x,y,o,xi,yi,r+, xi,zi,g* *,xi,wi,k.-),xi,wi,k.-) legend(legend(原始点原始点,线性点线性点,三次样条三次样条,三次多项式三次多项式)x0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.5221u二维数据插值二维数据插值 两个自变量两个自变量u在在MATLAB中，提供了解决二维插值问题的函数中，提供了解决二维插值问题的函数interp2，其调用格式为：，其调用格式为： Z1=interp2(X,

24、Y,Z,X1,Y1,method) X,Y是两个向量，分别描述两个参数的采样点是两个向量，分别描述两个参数的采样点 Z是与参数采样点对应的函数值是与参数采样点对应的函数值 X1,Y1是两个向量或标量，描述欲插值的点是两个向量或标量，描述欲插值的点 Z1是根据相应的插值方法得到的插值结果是根据相应的插值方法得到的插值结果 method的取值与一维插值函数相同的取值与一维插值函数相同22u例：已知某处山区地形选点测量坐标数据为：例：已知某处山区地形选点测量坐标数据为：x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

25、y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6u海拔高度数据为：海拔高度数据为：z=89 90 87 85 92 91 96 93 90 87 82z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 96 98 95 92 90 88 85 84 83 81 85 80 8

26、1 82 89 95 96 93 92 89 86 86 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 92 96 9

27、7 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 80 82 81 84 85 86 83 82 81 80 82 87 88 89 9

28、8 99 97 96 98 94 92 87 87 88 89 98 99 97 96 98 94 92 87u求对数据插值加密后形成的地貌图求对数据插值加密后形成的地貌图 23u绘原始数据图绘原始数据图x=0:.5:5;x=0:.5:5;y=0:.5:6;y=0:.5:6;z=89 90 87 85 92 91 96 93 90 87 82z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83

29、81 85 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 88 92 93 94 95 89 87 86 83 81

30、92 92 96 97 98 96 93 95 84 82 81 84 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82

31、80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87; 87 88 89 98 99 97 96 98 94 92 87;mesh(x,y,z)mesh(x,y,z)01234502468085909510024u绘插值后的地貌图绘插值后的地貌图 xi=linspace(0,5,50); xi=linspace(0,5,50); % %加密横坐标数据到加密横坐标数据到5050个个 yi=linspace(0,6,80); yi=linspace(0,6,80); % %加密纵坐标数据到加密纵坐标数据到6060个个 x

32、ii,yii=meshgrid(xixii,yii=meshgrid(xi,yi); %,yi); %生成网格数据生成网格数据 zii=interp2(x,y,z,xiizii=interp2(x,y,z,xii,yii,cubic); %,yii,cubic); %插值插值 mesh(xii,yii,zii) %mesh(xii,yii,zii) %加密后的地貌图加密后的地貌图0123450246758085909510025u比较比较 hold on % hold on % 保保持图形持图形 xx,yy=meshgridxx,yy=meshgrid(x,y); %(x,y); %生成网生成

33、网格数据格数据 plot3(xx,yy,z,oplot3(xx,yy,z,ob) b) % %原始数据用原始数据用OO绘出绘出 0123450246758085909510026u通常实验或测量得到的是一组离通常实验或测量得到的是一组离散数据序列（散数据序列（x xi i ,y,yi i）。）。u 如果数据序列（如果数据序列（x xi i ,y ,yi i）( (为一为一般起见般起见), ), i i=1,2, =1,2, ,m ,m ,含有不含有不可避免的误差（或称可避免的误差（或称“噪声噪声” ” ）u或者无法同时满足某特定的函数或者无法同时满足某特定的函数, ,那么，只能要求所作逼近函数

34、那么，只能要求所作逼近函数（x x）最优地靠近样点，即向量）最优地靠近样点，即向量Q=Q=（(x x1 1), (), (x x2 2) ), , , , ( (x xm)m)）T T与与Y=(Y=(y y1 1,y,y2 2, , ,y,ym m) )T T的误差或的误差或距离最小。距离最小。u按按Q Q与与Y Y之间误差最小原则作为之间误差最小原则作为“最优最优”标准构造的逼近函数，称标准构造的逼近函数，称为拟合函数。为拟合函数。 -202468101214161820050100150200Y X 024681005101520Y X 27u拟合与插值拟合与插值 共同点共同点都是通过已知

35、一些离散点集都是通过已知一些离散点集M M上的约束，求取一个定义在连续上的约束，求取一个定义在连续集合集合S(MS(M包含于包含于S)S)的未知连续函数，从而达到获取整体规律的的未知连续函数，从而达到获取整体规律的 目的，即通过目的，即通过“窥几斑窥几斑”来达到来达到“知全豹知全豹”。 不同点不同点从几何意义上讲，拟合是给定了空间中的一些点，找到一个从几何意义上讲，拟合是给定了空间中的一些点，找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点已知形式未知参数的连续曲面来最大限度地逼近这些点而插值是找到一个而插值是找到一个( ( 或几个分片光滑的或几个分片光滑的) )连续曲面来穿过这些连续

36、曲面来穿过这些点。点。 28u最常用的的曲线拟合是最常用的的曲线拟合是最小二乘法曲线拟合最小二乘法曲线拟合，拟合结果可，拟合结果可使误差的平方和最小。使误差的平方和最小。u在在MATLAB中，用中，用polyfit函数来求得最小二乘拟合多项式函数来求得最小二乘拟合多项式的系数，再用的系数，再用polyval函数按所得的多项式计算所给出的点函数按所得的多项式计算所给出的点上的函数近似值。上的函数近似值。 polyfit函数的调用格式为：函数的调用格式为： P,S=polyfit(X,Y,m)函数根据采样点函数根据采样点X和采样点函数值和采样点函数值Y，产生一个，产生一个m次多项式次多项式P及其在

37、及其在采样点的误差向量采样点的误差向量S。其中其中X,Y是两个等长的向量，是两个等长的向量，P是一个长度为是一个长度为m+1的向量，的向量，P的元素的元素为多项式系数。为多项式系数。 polyval函数的功能是按多项式的系数计算函数的功能是按多项式的系数计算x点多项式的值。点多项式的值。29u例例 用一个用一个3次多项式在区间次多项式在区间0，2内逼近函数。内逼近函数。 X=linspace(0,2*pi,50); Y=sin(X); P=polyfit(X,Y,3) %得到得到3次多项式的系数和误差次多项式的系数和误差 X=linspace(0,2X=linspace(0,2* *pi,20

38、);pi,20); Y=sin(X);Y=sin(X); Y1=polyval(P,X)Y1=polyval(P,X) plot(X,Y,:o,X,Y1,-plot(X,Y,:o,X,Y1,-* *)01234567-1-0.8-0.6-0.4-0.200.20.40.60.8130uN N次多项式表示为次多项式表示为 P(x)=aP(x)=a0 0 x xn n+a+a1 1x xn-1n-1+a+a2 2x xn-2n-2a an-1n-1x+ax+an nuMatlabMatlab中中n n次多项式用一个长度为次多项式用一个长度为n+1n+1的行向量（系数向量的行向量（系数向量）表示）表

39、示 aa0 0 a a1 1 a an-1n-1 a an n u多项式的四则运算多项式的四则运算 多项式的加减运算多项式的加减运算系数向量的加减运算系数向量的加减运算要求次数相同，不足时用要求次数相同，不足时用“0”补起补起例例 54322( )352756( )353f xxxxxxg xxx31u多项式乘法运算多项式乘法运算 函数函数conv(P1,P2)用于求多项式用于求多项式P1和和P2的乘积。这里，的乘积。这里，P1、P2是两个多项式系数向是两个多项式系数向量量u多项式除法运算多项式除法运算 函数函数Q,r=deconv(P1,P2)用于对多项式用于对多项式P1和和P2作除法运算。

40、其中作除法运算。其中Q返回多项式返回多项式P1除除以以P2的商式，的商式，r返回返回P1除以除以P2的余式。这里，的余式。这里，Q和和r仍是多项式系数向量。仍是多项式系数向量。 deconv是是conv的逆函数，即有的逆函数，即有P1=conv(P2,Q)+r。u例例 求求f(x)+g(x)f(x)+g(x)、f(x)-g(x)f(x)-g(x)。 求求f(x)f(x)g(x)g(x)、f(x)/g(x)f(x)/g(x)。 f=3,-5,2,-7,5,6;g=3,5,-3;g1=0,0,0,g;f=3,-5,2,-7,5,6;g=3,5,-3;g1=0,0,0,g; f+g1 %f+g1 %

41、求求f(x)+g(x)f(x)+g(x) f-g1 %f-g1 %求求f(x)-g(x)f(x)-g(x) conv(f,g) %conv(f,g) %求求f(x)f(x)* *g(x)g(x) Q,r=deconv(f,g) %Q,r=deconv(f,g) %求求f(x)/g(x)f(x)/g(x)，商式送，商式送Q Q，余式送，余式送r r。54322( )352756( )353f xxxxxxg xxx32u多项式的导函数多项式的导函数 p=polyder(P)：求多项式：求多项式P的导函数的导函数 p=polyder(P,Q)：求：求PQ的导函数的导函数 p,q=polyder(P

42、,Q)：求：求P/Q的导函数，导函数的分子存入的导函数，导函数的分子存入p，分母存入，分母存入qu例例 求有理分式的导数。求有理分式的导数。 P=3,5,0,-8,1,-5;P=3,5,0,-8,1,-5; Q=10,5,0,0,6,0,0,7,-1,0,-100;Q=10,5,0,0,6,0,0,7,-1,0,-100; p,q=polyder(P,Q)p,q=polyder(P,Q)5421096323585( )10567100 xxxxf xxxxxx33u多项式求值多项式求值MATLAB提供了两种求多项式值的函数：提供了两种求多项式值的函数：polyval与与polyvalm，它们的

43、输入参数均，它们的输入参数均为多项式系数向量为多项式系数向量P和自变量和自变量x。两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。u代数多项式求值代数多项式求值Y=polyval(P,x)若若x为一数值，则求多项式在该点的值；若为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。每个元素求其多项式的值。u矩阵多项式求值矩阵多项式求值polyvalm函数用来求矩阵多项式的值，其调用格式与函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同

44、。相同，但含义不同。polyvalm函数要求函数要求x为方阵，它以方阵为自变量求多项式的值。为方阵，它以方阵为自变量求多项式的值。upolyval与与polyvalm的区别的区别设设A为方阵，为方阵，P代表多项式代表多项式x3-5x2+8，那么，那么polyvalm(P,A)的含义是：的含义是：A*A*A-5*A*A+8*eye(size(A)polyval(P,A)的含义是：的含义是：A.*A.*A-5*A.*A+8*ones(size(A)34u以多项式以多项式x4+8x3-10为例，取一个为例，取一个22矩阵为自变量分别用矩阵为自变量分别用polyval和和polyvalm计算该多项式的

45、值。计算该多项式的值。 A=1,8,0,0,-10; % A=1,8,0,0,-10; % 多项式系数多项式系数 x=-1,1.2;2,-1.8 % x=-1,1.2;2,-1.8 % 给出一个矩阵给出一个矩阵x x y1=polyval(A,x) % y1=polyval(A,x) % 计算代数多项式的值计算代数多项式的值 y2=polyvalm(A,x) % y2=polyvalm(A,x) % 计算矩阵多项式的值计算矩阵多项式的值35u多项式求根多项式求根un次多项式具有次多项式具有n个根，当然这些根可能是实根，也可能含个根，当然这些根可能是实根，也可能含有若干对共轭复根有若干对共轭复根

46、 MATLAB提供的提供的roots函数用于求多项式的全部根，其调用格式为函数用于求多项式的全部根，其调用格式为：x=roots(P)其中其中P为多项式的系数向量，求得的根赋给向量为多项式的系数向量，求得的根赋给向量x，即，即x(1),x(2),x(n)分别代表多项式的分别代表多项式的n个根。个根。 P=poly(x)若若x为具有为具有n个元素的向量，则个元素的向量，则poly(x)建立以建立以x为其根的多项式，且将为其根的多项式，且将该多项式的系数赋给向量该多项式的系数赋给向量P。36u例例 已知已知 (1) 计算计算f(x)=0 的全部根。的全部根。(2) 由方程由方程f(x)=0的根构造

47、一个多项式的根构造一个多项式g(x)，并与，并与f(x)进行对进行对比。比。命令如下：命令如下：P=3,0,4,-5,-7.2,5;X=roots(P) %求方程求方程f(x)=0的根的根G=poly(X) %求多项式求多项式g(x)52 . 7543)(235xxxxxf37u数据处理与多项式计算数据处理与多项式计算u数值微积分数值微积分u离散傅立叶变换离散傅立叶变换u线性方程组求解线性方程组求解u非线性方程与最优化问题求解非线性方程与最优化问题求解u常微分方程的数值求解常微分方程的数值求解u稀疏矩阵稀疏矩阵38u数值微分的两种方法数值微分的两种方法用多项式或样条函数用多项式或样条函数g(x

48、)对对f(x)进行逼近（插值或拟合），然后用逼近函数进行逼近（插值或拟合），然后用逼近函数g(x)在点在点x处的导数作为处的导数作为f(x)在点在点x处的导数处的导数用用f(x)在点在点x处的某种差商作为其导数处的某种差商作为其导数u在在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其，其调用格式为：调用格式为：DX=diff(X)：计算向量：计算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1DX=diff(X,n)：计算：计算X的的n阶向前差分。例如，阶向前差分。例如

49、，diff(X,2)=diff(diff(X)DX=diff(A,n,dim)：计算矩阵：计算矩阵A的的n阶差分，阶差分，dim=1时时(缺省状态缺省状态)，按列计算差分；，按列计算差分；dim=2，按行计算差分。按行计算差分。u例例 设设x由由0,2间均匀分布的间均匀分布的10个点组成，求个点组成，求sinx的的13阶差分阶差分X=linspace(0,2*pi,10);Y=sin(X);DY=diff(Y); %计算计算Y的一阶差分的一阶差分D2Y=diff(Y,2); %计算计算Y的二阶差分，也可用命令的二阶差分，也可用命令diff(DY)计算计算D3Y=diff(Y,3); %计算计算

50、Y的三阶差分，也可用的三阶差分，也可用diff(D2Y)或或diff(DY,2)39u数值积分基本原理数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson) 法、牛法、牛顿柯特斯顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整法等都是经常采用的方法。它们的基本思想都是将整个积分区间个积分区间a,b分成分成n个子区间个子区间xi,xi+1，i=1,2,n，其中，其中x1=a，xn+1=b。这样求定。这样求定积分问题就分解为求和问题。积分问题就分解为求和问题。u数值积分的实现数值积

51、分的实现变步长辛普生法变步长辛普生法 I,n=quad(fname,a,b,tol,trace) I,n=quad(fname,a,b,tol,trace) 其中其中fnamefname是被积函数名。是被积函数名。a a和和b b分别是定积分的下限和上限。分别是定积分的下限和上限。toltol用来控制积分精用来控制积分精度，缺省时取度，缺省时取tol=0.001tol=0.001。tracetrace控制是否展现积分过程，若取非控制是否展现积分过程，若取非0 0则展现积分过则展现积分过程，取程，取0 0则不展现，缺省时取则不展现，缺省时取trace=0trace=0。返回参数。返回参数I I即

52、定积分值，即定积分值，n n为被积函数的调为被积函数的调用次数用次数牛顿牛顿- -柯斯特法柯斯特法I,n=quad8(fname,a,b,tol,trace)I,n=quad8(fname,a,b,tol,trace)LobbatoLobbato法法高精度高精度I,n=quadl(fname,a,b,tol,trace)I,n=quadl(fname,a,b,tol,trace)被积函数由一个表格定义被积函数由一个表格定义trapz(X,Y) trapz(X,Y) 其中其中向量向量X、Y 分别代表分别代表样本点和样本值，定义函数关系样本点和样本值，定义函数关系Y=f(X) 二重定积分二重定积分

53、I=dblquad(f,a,b,c,d,tol,trace)I=dblquad(f,a,b,c,d,tol,trace)40uFnameFname为描述被为描述被积函数的字符串变量，定义有两种方式积函数的字符串变量，定义有两种方式 建立关于被积函数的函数文件建立关于被积函数的函数文件 function ex=ex(x) ex=exp(-x.2); 使用语句函数使用语句函数(内联函数内联函数) Inline() Inline()，第一个变量为，第一个变量为被被积函数本身，第二个变量为自变量，也积函数本身，第二个变量为自变量，也可以带多个自变量可以带多个自变量 g=inline(exp(-x.2)

54、; %定义一个语句函数定义一个语句函数g(x)=exp(-x2)41u例例 用两种不同的方法求定积分用两种不同的方法求定积分方法一：先建立一个函数文件方法一：先建立一个函数文件ex.m：function ex=ex(x)ex=exp(-x.2);然后在然后在MATLAB命令窗口，输入命令：命令窗口，输入命令：format longI=quad(ex,0,1) %注意函数名应加字符引号注意函数名应加字符引号I = 0.74682418072642I=quadl(ex,0,1)I = 0.74682413398845方法二：也可不建立关于被积函数的函数文件，而使用语句函数方法二：也可不建立关于被积

55、函数的函数文件，而使用语句函数(内联函数内联函数)求解求解，命令如下：，命令如下：g=inline(exp(-x.2); %定义一个语句函数定义一个语句函数g(x)=exp(-x2)I=quadl(g,0,1) %注意函数名不加注意函数名不加号号I = 0.74682413398845format short102dxeIx42u数据处理与多项式计算数据处理与多项式计算u数值微积分数值微积分u离散傅立叶变换离散傅立叶变换u线性方程组求解线性方程组求解u非线性方程与最优化问题求解非线性方程与最优化问题求解u常微分方程的数值求解常微分方程的数值求解u稀疏矩阵稀疏矩阵43u傅立叶变换傅立叶变换信号信

56、号以时间为自变量以频率为自变量傅里叶变换逆傅里叶变换44u离散傅立叶变换算法离散傅立叶变换算法 在计算机上实现信在计算机上实现信号的频谱分析要求号的频谱分析要求时域和频域都是时域和频域都是离散的；离散的；时域和频域都是时域和频域都是有限长的。有限长的。 有卷积波纹 N N 原函数 用于抽样 抽样后 用于截断 截断后 用于延拓 延拓后 定义DFT 叠加干涉 0 t 0 t 0 t 0 t 0 t 0 t 0 t 0 t 或nTs 0 f或kf0 0 f 0 f 0 f 0 f 0 f 0 f 0 f 45u离散傅立叶变换的实现离散傅立叶变换的实现采用快速采用快速傅立叶变换（傅立叶变换（FFTFF

57、T）FFTFFT不是一种新的变换，而是不是一种新的变换，而是DFTDFT的快速算法。的快速算法。u一维离散傅立叶变换函数，其调用格式与功能为：一维离散傅立叶变换函数，其调用格式与功能为：(1)fft(X)：返回向量：返回向量X的离散傅立叶变换。设的离散傅立叶变换。设X的长度的长度(即元素个数即元素个数)为为N，若，若N为为2的幂次，则为以的幂次，则为以2为基数的快速傅立叶变换，否则为运算速度很慢的非为基数的快速傅立叶变换，否则为运算速度很慢的非2幂次的算法。对于矩阵幂次的算法。对于矩阵X，fft(X)应用于矩阵的每一列。应用于矩阵的每一列。(2) fft(X,N)：计算：计算N点离散傅立叶变换

58、。它限定向量的长度为点离散傅立叶变换。它限定向量的长度为N，若，若X的长度小的长度小于于N，则不足部分补上零；若大于，则不足部分补上零；若大于N，则删去超出，则删去超出N的那些元素。对于矩阵的那些元素。对于矩阵X，它同样应用于矩阵的每一列，只是限定了向量的长度为，它同样应用于矩阵的每一列，只是限定了向量的长度为N。(3) fft(X,dim)或或fft(X,N,dim)：这是对于矩阵而言的函数调用格式，前者的功：这是对于矩阵而言的函数调用格式，前者的功能与能与FFT(X)基本相同，而后者则与基本相同，而后者则与FFT(X,N)基本相同。只是当参数基本相同。只是当参数dim=1时，该函数作用于时

59、，该函数作用于X的每一列；当的每一列；当dim=2时，则作用于时，则作用于X的每一行。的每一行。46u一维离散傅立叶逆变换函数一维离散傅立叶逆变换函数 一维离散傅立叶逆变换函数是一维离散傅立叶逆变换函数是ifft。 ifft(F)返回返回F的一维离散傅立叶逆变换；的一维离散傅立叶逆变换； ifft(F,N)为为N点逆变换；点逆变换；ifft(F,dim)或或ifft(F,N,dim)则由则由N或或dim确确定逆变换的点数或操作方向。定逆变换的点数或操作方向。u例例 给定数学函数给定数学函数 x(t)=12sin(210t+/4)+5cos(240t) 取取N=128，试对，试对t从从01秒采样

60、，用秒采样，用fft作快速傅立叶变换，绘制相应作快速傅立叶变换，绘制相应的的 振幅振幅-频率图。频率图。 47u程序如下：程序如下：N=128; % 采样点数采样点数T=1; % 采样时间终点采样时间终点t=linspace(0,T,N); % 给出给出N个采样时间个采样时间ti(I=1:N)x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t); % 求各采样点样本值求各采样点样本值xdt=t(2)-t(1); % 采样周期采样周期f=1/dt; % 采样频率采样频率(Hz)X=fft(x); % 计算计算x的快速傅立叶变换的快速傅立叶变换XF=X(1:N/2+1