《第三章数值运算基础优秀PPT.ppt》由会员分享，可在线阅读，更多相关《第三章数值运算基础优秀PPT.ppt（40页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第三章 数值运算基础1现在学习的是第1页，共40页 （3 3）由根矢量创建多项式）由根矢量创建多项式 已知一个多项式的全部根已知一个多项式的全部根X X求多项式系数的函数是求多项式系数的函数是poly(X)poly(X)，该函数返回以该函数返回以X X为全部根的一个多项式为全部根的一个多项式P P，当，当X X是一个长度为是一个长度为m m的向量时，的向量时，P P是一个长度为是一个长度为m+1m+1的向量。的向量。（2）特征多项式输入法）特征多项式输入法n阶方阵的特征多项式系数矢量一定是阶方阵的特征多项式系数矢量一定是n+1阶的，同时特征阶的，同时特征多项式系数矢量的第一个元素必须为多项式系
2、数矢量的第一个元素必须为1。例、求矩阵例、求矩阵A=123；456；789的特征多项式系数，的特征多项式系数，并转换为多项式形式。并转换为多项式形式。A=123；456；789；P=poly（A）f=poly2sym（P）例、由根矢量例、由根矢量-0.5-0.3+0.4i-0.3-0.4i创建多项式创建多项式R=-0.5-0.3+0.4i-0.3-0.4i;P=poly（R）f=poly2sym（P）现在学习的是第2页，共40页二、多项式运算二、多项式运算1、求多项式的值、求多项式的值求多项式求多项式p(x)在某点或某些点的函数值的函数是在某点或某些点的函数值的函数是polyval(P,x)。
3、若若x为一数值，则求多项式在该点的值；为一数值，则求多项式在该点的值；若若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值;若若x为方阵，则可用函数为方阵，则可用函数polyvalm（P，x）按矩阵的运算规则求多项式值。按矩阵的运算规则求多项式值。例、（例、（1）求多项式）求多项式3x2+2x+1在在2、5、7、9处的值；处的值；（2）求多项式）求多项式3x2+2x+1对于矩阵对于矩阵25；79的值。的值。P=321；polyval（P，2579）polyvalm（P，25；79）polyval（P，25；79）现在学习的是第3页，
4、共40页2、求多项式的根、求多项式的根求多项式求多项式p(x)的根的函数是的根的函数是roots(P)，这里，这里，P是是p(x)的的系数向量，该函数返回方程系数向量，该函数返回方程p(x)=0的全部根的全部根(含重根，复根含重根，复根)。例、求多项式例、求多项式x3-27的根的根A=100-27；r=roots（A）在在MATLAB中约定：中约定：3 3、多项式的四则运算、多项式的四则运算(1)(1)多项式的加减法：对应系数矢量相加减多项式的加减法：对应系数矢量相加减(2)(2)多项式的乘法多项式的乘法 函数函数conv(P1,P2)conv(P1,P2)用于求多项式用于求多项式P1P1和和
5、P2P2的乘积；的乘积；也是矢量的卷积函数也是矢量的卷积函数现在学习的是第4页，共40页(3)多项式的除法多项式的除法函数函数Q,r=deconv(P1,P2)用于对多项式用于对多项式P1和和P2作除法作除法运算。其中运算。其中Q返回多项式返回多项式P1除以除以P2的商式，的商式，r返回返回P1除以除以P2的余式。的余式。这里，这里，Q和和r仍是多项式系数向量。仍是多项式系数向量。deconv是是conv的逆函数，即有的逆函数，即有P1=conv(P2,Q)+r。例例1、求多项式、求多项式x3+3x2-2x+2和和x3+2x2+x-3的乘积。的乘积。P1=13-22；P2=121-3；P=co
6、nv（P1，P2）例例2、现在学习的是第5页，共40页f1=121；k1=-101；f2=ones（1，5）；）；k2=-2：2f=conv（f1，f2）4、多项式的微积分、多项式的微积分多项式的微分由函数多项式的微分由函数polyder实现实现多项式的积分由函数多项式的积分由函数polyint实现实现例、计算多项式例、计算多项式5x6+3x4-12x3+1的微分的微分P=503-12001；A=polyder（P）poly2sym（A）B=polyint（A）poly2sym（B）5、多项式的部分分式展开、多项式的部分分式展开（1）对于多项式）对于多项式b（x）和不含重根的）和不含重根的n阶
7、多项式阶多项式a（x）之比，有如下展开：）之比，有如下展开：现在学习的是第6页，共40页（2）当）当a（x）有）有m个重根个重根p（j）=p（j+1）=p（j+m-1）时，）时，相应的部分分式为：相应的部分分式为：在在MATLAB中，利用函数中，利用函数residue，调用格式为调用格式为r，p，k=residue（b，a）留数留数极点极点直项直项分子分子分母分母逆函数：逆函数：b，a=residue（r，p，k）现在学习的是第7页，共40页例例1、求、求部分分式展开。部分分式展开。A=-4083；B=53-27；r，p，k=residue（b，a）例例2、求解微分方程：、求解微分方程：现在学
8、习的是第8页，共40页MATLAB实现：实现：B=1，2；A=1，4，3；R，P，K=residue（B，A）得到：得到：R=0.5P=-3K=0.5-1即有：即有：MATLAB实现：实现：A1=conv（1，4，3，1，1）；R1，P1，K1=residue（B，A1）现在学习的是第9页，共40页得到：得到：R1=-0.25P1=-3K1=0.25-10.5-1即有：即有：现在学习的是第10页，共40页6、多项式拟合、多项式拟合什么是多项式拟合？什么是多项式拟合？多项式拟合函数：多项式拟合函数：p=ployfit(x,y,n)应用最小二乘法应用最小二乘法求出求出n阶拟合多项式阶拟合多项式p(
9、x)，即用，即用p(x(i)拟合拟合y(i).p,S=ployfit(x,y,n)返回返回n阶拟合多项式阶拟合多项式p(x)合包合包括误差估计的结构括误差估计的结构S。p,S，mu=ployfit(x,y,n)在这种用法中，用在这种用法中，用代替代替x，求拟合多项式，求拟合多项式p11现在学习的是第11页，共40页3.2线性代数线性代数在在分分析析及及解解决决问问题题的的过过程程中中，通通常常根根据据已已知知条条件件尝尝试试将将相相同同以以方方程程的的形形式式来来表表示示，再再由由求求出出的的方方程程解解来来进进一一步步了了解解相相同，所以解方程是非常重要的。同，所以解方程是非常重要的。解解线
10、线性性方方程程就就是是找找出出是是否否存存在在一一个个唯唯一一矩阵矩阵x,使得矩阵使得矩阵a、b满足一些关系：满足一些关系：ax=b或或xa=b12现在学习的是第12页，共40页MATLAB以斜线和反斜线来表示除法运算因子，其中：以斜线和反斜线来表示除法运算因子，其中：x=ab=a-1*b是方程式是方程式ax=b的解。的解。x=b/a=b*a-1是方程式是方程式xa=b的解。的解。可能由以下可能由以下3种情况：种情况：（1）m=n（方阵系统），（方阵系统），可以尝试计算精确解。可以尝试计算精确解。（2）mn（超定系统），可以尝试计算最小二乘解。（超定系统），可以尝试计算最小二乘解。mn（欠定系
11、统），（欠定系统），可以尝试计算含有最小可以尝试计算含有最小m的解。的解。13现在学习的是第13页，共40页反反斜斜线线运运算算因因子子第第一一不不同同形形式式的的参参数数矩矩阵阵，采采 用用 不不 同同 的的 运运 算算 法法 则则 来来 处处 理理。Matlab会会自自动动检检测测参参数数矩矩阵阵，以以区区分分下下面几种形式：面几种形式：三角矩阵三角矩阵对称正定矩阵对称正定矩阵非奇异矩阵非奇异矩阵超定系统超定系统欠定系统欠定系统14现在学习的是第14页，共40页一、一、方阵系统方阵系统最最常常见见的的线线性性方方程程是是系系数数矩矩阵阵为为方方阵阵a和和由由常常数数项项组组成成列列矢矢量量
12、b的的情情况况，则则x可可写写成成x=ab，其中，其中x和和b的尺寸相同。的尺寸相同。【例】【例】求方阵系统的根。求方阵系统的根。在在命命令令窗窗口口产产生生随随机机方方阵阵a和和矢矢量量b，计计算算ab:a=fix(15*rand(3,3)b=fix(15*rand(3,1)x=ab15现在学习的是第15页，共40页【例】【例】a和和b均为方阵，求方阵系统的根。均为方阵，求方阵系统的根。在命令窗口产生随机方阵在命令窗口产生随机方阵a和和b，计算，计算ab：a=fix(15*rand(3,3)b=fix(15*rand(3,3)x=ab假假如如方方阵阵a的的各各个个行行矢矢量量线线性性相相关关
13、，则则方方阵阵a为为奇奇异阵，这也使线性方阵异阵，这也使线性方阵ax=b有无穷多组解。有无穷多组解。假假如如方方阵阵a近近似似奇奇异异矩矩阵阵，则则反反斜斜线线运运算算因因子子将将发发出出警警告告信信息息。假假如如方方阵阵a确确定定是是奇奇异异矩矩阵阵，则则反反斜斜线线运运算算因子将发出警告信息。因子将发出警告信息。16现在学习的是第16页，共40页3.3数据分析数据分析MATLAB对数据分析命令有两条约定：对数据分析命令有两条约定：（1）若输入量）若输入量X为矢量，则不论是行矢量还是列矢量，运算是为矢量，则不论是行矢量还是列矢量，运算是对整个矢量进行，其结果为一个数值；对整个矢量进行，其结果
14、为一个数值；（2）若输入量）若输入量X为矩阵，则函数运算按列进行，其结果为一个为矩阵，则函数运算按列进行，其结果为一个行矢量。行矢量。一、数据统计与分析一、数据统计与分析1.求矩阵最大和最小元素求矩阵最大和最小元素(1)求向量的最大元素求向量的最大元素y=max(X)返回向量返回向量X的最大元素存入的最大元素存入y。y,I=max(X)返回向量返回向量X的最大元素存入的最大元素存入y，最大元素的序号存入，最大元素的序号存入I。(2)求矩阵的最大元素求矩阵的最大元素max(A)返回一个行向量，向量的第返回一个行向量，向量的第i个元素是个元素是A矩阵的第矩阵的第i列上列上的最大元素。的最大元素。现
15、在学习的是第17页，共40页Y,U=max(A)返回两个行向量，返回两个行向量，Y向量记录向量记录A的每列的最大元素，的每列的最大元素，U向量记录每列最大元素的行号。向量记录每列最大元素的行号。max(A,dim)dim取取1或或2。dim取取1时，该函数和时，该函数和max(A)完全完全相同。相同。dim取取2时，该函数返回一个列向量，其第时，该函数返回一个列向量，其第i个元素是个元素是A矩阵矩阵的第的第i行上的最大元素。行上的最大元素。（3）求最小值函数：）求最小值函数：min（）用法同用法同max（）例、例、求矩阵求矩阵A的每行及每列的最大和最小元素，并求整个矩阵的最大和最小的每行及每列
16、的最大和最小元素，并求整个矩阵的最大和最小元。元。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)%求整个矩阵的最小元素求整个矩阵的最小元素现在学习的是第18页，共40页2.求矩阵的平均值和中值求矩阵的平均值和中值求矩阵和向量元素的平均值的函数是求矩阵和向量元素的平均值的函数是mean，求矩阵和向量元素的中值
17、的函数是求矩阵和向量元素的中值的函数是median。其调用方法：其调用方法：y=mean(X)返回向量返回向量X的平均值存入的平均值存入y。mean(A)返回一个行向量，向量的第返回一个行向量，向量的第i个个元素是元素是A矩阵的第矩阵的第i列上的元素的平均值。列上的元素的平均值。mean(A，2)返回一个列向量，向量的第返回一个列向量，向量的第i个个元素是元素是A矩阵的第行上的元素的平均值。矩阵的第行上的元素的平均值。3.矩阵元素求和与求积矩阵元素求和与求积矩阵和向量求和的基本函数是矩阵和向量求和的基本函数是sum矩阵和向量求积的基本函数是矩阵和向量求积的基本函数是prod，其使用方法与其使用
18、方法与mean（）函数类似）函数类似例、求矩阵例、求矩阵A的每行元素的乘积和全部元素的乘积的每行元素的乘积和全部元素的乘积A=1,2,3,4;5,6,7,8;9,10,11,12;S=prod(A,2)prod(S)%求求A的全部元素的乘积的全部元素的乘积现在学习的是第19页，共40页4.矩阵元素累加和与累乘积矩阵元素累加和与累乘积求向量和矩阵元素的累加和求向量和矩阵元素的累加和cumsum求向量和矩阵元素的累乘积求向量和矩阵元素的累乘积cumprod其使用方法与其使用方法与mean（）函数类似）函数类似例、求向量例、求向量X=(1！，！，2！，！，3！，！，10！)X=cumprod(1:1
19、0)5.标准方差标准方差MATLAB中，提供了计算数据序列的标准方差的函数中，提供了计算数据序列的标准方差的函数std。对于向量对于向量X，std(X)返回一个标准方差。返回一个标准方差。对于矩阵对于矩阵A，std(A)返回一个行向量，是矩阵返回一个行向量，是矩阵A各列的标准方差。各列的标准方差。A=12345；678910；1112131415；B=std（A）现在学习的是第20页，共40页6.元素递增排序元素递增排序对向量对向量X排序函数是排序函数是sort(X)，返回一个对，返回一个对X中的元素按升序排列的新中的元素按升序排列的新向量。向量。对矩阵对矩阵A的各列的各列(或行或行)重新排序
20、，其调用格式重新排序，其调用格式Y,I=sort(A,dim)其中其中dim指明对指明对A的列还是行进行排序，若的列还是行进行排序，若dim=1，则按列排，若，则按列排，若dim=2，则按行排。，则按行排。Y是排序后的矩阵，而是排序后的矩阵，而I记录记录Y中的元素在中的元素在A中位置。中位置。例、对矩阵做各种排序。例、对矩阵做各种排序。A=1,-8,5;4,12,6;13,7,-13;sort(A)%对对A的每列按升序排序的每列按升序排序sort(A,2)%对对A的每行按升序排序的每行按升序排序X,I=sort(A)%对对A按列排序，并将每个元素所在行号送矩阵按列排序，并将每个元素所在行号送矩
21、阵I现在学习的是第21页，共40页二、协方差阵二、协方差阵(Covariancematrix)和和相关系数阵相关系数阵(Correlationcoefficients)1、函数、函数cov（x，y）求两随机变量的协方差矩阵求两随机变量的协方差矩阵2、corrcoef（x，y）求两随机变量的相关系数矩阵求两随机变量的相关系数矩阵现在学习的是第22页，共40页三、插值三、插值插值（插值（interpolation）是在已知数据之间计算估计值的过程，在数字信号）是在已知数据之间计算估计值的过程，在数字信号处理中应用广泛。处理中应用广泛。1.一维数值插值一维数值插值（1）interp1函数调用格式为：
22、函数调用格式为：Y1=interp1(X,Y,X1,method)函数根据函数根据X、Y的值，计算函数在的值，计算函数在X1处的值。处的值。X、Y是两个等长的已知向量，分别描述采样点和样本值，是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，是一个向量或标量，描述欲插值的点，Y1是一个与是一个与X1等长的插值结果。等长的插值结果。（2）method是插值方法，允许的取值有是插值方法，允许的取值有linear(线性插值线性插值):在两采样点间连直线，插值在两采样点间连直线，插值X1对应的对应的Y1在直线上在直线上nearest(最近插值最近插值)：将插值结果设置
23、为最近数据点的值；：将插值结果设置为最近数据点的值；spline(三次样条插值三次样条插值)：通过数据点拟合出三次样条曲线：通过数据点拟合出三次样条曲线cubic（三次多项式插值）：通过分段立方插值法计算插值结果；（三次多项式插值）：通过分段立方插值法计算插值结果；现在学习的是第23页，共40页例、用不同的插值方法计算例、用不同的插值方法计算sin(x)在在/2点的值。点的值。这是一个一维插值问题。这是一个一维插值问题。X=0:0.2:pi;Y=sin(X);%给出给出X、Yinterp1(X,Y,pi/2)%用缺省方法用缺省方法(即线性插值法即线性插值法)计算计算sin(/2)interp1
24、(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)（3）MATLAB中有一个专门的三次样条插值函中有一个专门的三次样条插值函Y1=spline(X,Y,X1)，其功能及使用方法与函数其功能及使用方法与函数Y1=interp1(X,Y,X1,spline)完全相同。完全相同。例
25、、例、X=0：10；Y=sin（X）；）；X1=0：0.5：10；Y1=spline（X，Y，X1）；）；plot（X，Y，k，X1，Y1，*r）现在学习的是第24页，共40页2.二维数值插值二维数值插值Z1=interp2(X,Y,Z,X1,Y1,method)3.三维数值插值三维数值插值W1=interp3(X,Y,Z,W,X1,Y1,Z1,method)4、选择插值法时考虑的因素：、选择插值法时考虑的因素：运算时间、占用内存、曲线的光滑程度运算时间、占用内存、曲线的光滑程度5、信号的整数倍插值、信号的整数倍插值信号插值，是一个增大采样频率来增加数据的过程。然而信号插值，是一个增大采样频率
26、来增加数据的过程。然而这样的简单插值会在频域内产生一些多余图像，所以在插值这样的简单插值会在频域内产生一些多余图像，所以在插值后要用低通滤波器去掉多余的频域图像。后要用低通滤波器去掉多余的频域图像。函数调用方式为：函数调用方式为：现在学习的是第25页，共40页（1）y=interp（x，r）：将信号）：将信号x的采样频率提高为原来的的采样频率提高为原来的r倍。倍。在插值后，默认地采用了对称型的截频低通滤波器。在插值后，默认地采用了对称型的截频低通滤波器。（2）y=interp（x，r，l，alpha）：参数）：参数l和参数和参数alpha用于指定所用于指定所采用的截频低通滤波器的长度和截止频率
27、，默认值分别为采用的截频低通滤波器的长度和截止频率，默认值分别为4和和0.5。（3）y，b=interp（x，r，l，alpha）：）：b为截频低通滤波器的系为截频低通滤波器的系数向量的返回值。数向量的返回值。例、对原信号进行例、对原信号进行4倍插值。倍插值。6、信号的整数倍抽取、信号的整数倍抽取信号抽取（信号抽取（decimation），是一个减小采样频率来去掉多余数据），是一个减小采样频率来去掉多余数据的过程。然而这样的抽样不满组香农定理时，抽样信号会发生混叠。的过程。然而这样的抽样不满组香农定理时，抽样信号会发生混叠。所以在抽样前，必须对被采样信号做低通滤波来压缩频带，然后再所以在抽样前
28、，必须对被采样信号做低通滤波来压缩频带，然后再抽取，以避免混叠现象。抽取，以避免混叠现象。现在学习的是第26页，共40页（3）y=decimate（x，r，fir）：采用）：采用30点的点的FIR型滤波器来型滤波器来压缩频带。压缩频带。（4）y=decimate（x，r，n，fir）：采用）：采用n点的点的FIR型滤波器型滤波器来来压缩频带。压缩频带。例、对原信号进行例、对原信号进行4倍抽取。倍抽取。函数调用方式为：函数调用方式为：（1）y=decimate（x，r）：将信号）：将信号x的采样频率降低为原来的的采样频率降低为原来的1/r在抽取之前，默认地采用了在抽取之前，默认地采用了8阶阶Ch
29、eyshev型低通滤波器。型低通滤波器。（2）y=decimate（x，r，n）：参数）：参数n为采用的为采用的Cheyshev型低通型低通滤波器的阶数。滤波器的阶数。现在学习的是第27页，共40页*3.5*3.5 数字信号处理初步数字信号处理初步一、波形的产生1、周期波形（1）方波：square函数 X=square（t）：产生周期为2、峰值为1的方波X=square（t，duty）：产生指定周期、峰值为1的方波；duty 为占空比，即信号为正的周期所占的比例。例、产生周期为产生周期为1的标准方波和占空比为的标准方波和占空比为80%的方波。的方波。（2）锯齿波或三角波 sawtooth函数
30、X=sawtooth（t）：产生周期为2、峰值为1的锯齿波X=sawtooth（t，width）：根据width值产生不同形状的三角波 注：width为01之间的标量，指定一个周期内X最大值出现 的位置，width是位置横坐标与周期的比值。例、产生周期为产生周期为1的锯齿波波和对称三角波。的锯齿波波和对称三角波。现在学习的是第28页，共40页2、特殊波形（1）抽样函数波形 sinc函数X=sinc（t）：产生以t为变量的抽样函数波形例、t=-10：0.01：10；y=sinc（t）；plot（t，y）（2）单位脉冲序列使用MATLAB语言描述为：X=-7：7；Y=zeros（1，7），1，ze
31、ros（1，7）；stem（X，Y）axis（-5，5，0，2）现在学习的是第29页，共40页（3）单位阶跃序列使用MATLAB语言描述为：X=-7：7；Y=zeros（1，7），ones（1，8）；stem（X，Y，r*）axis（-7，7，0，2）（4）实指数序列使用MATLAB语言描述为：X=-7：7；Y=2.X；stem（X，Y，r*）axis（-7，7，0，50）现在学习的是第30页，共40页二、离散傅里叶变换 离散傅里叶变换是信号分析中的一种重要工具，它将时域内的问题转化为频域内的问题，在很多情况下大大地简化了问题的求解过程；另外计算时域信号的频谱也主要依靠离散傅里叶变换来完成。离
32、散序列x（k）的离散傅里叶变换（DFT）定义为：离散序列x（n）的离散傅里叶逆变换定义为：由于MATLAB中序列的下标不允许出现0，而是从1开始，因此在MATLAB中相应的傅立叶变换和逆变换的表达式为：现在学习的是第31页，共40页（1）fft函数：一维快速傅立叶变换 调用形式：Y=fft（X，N）返回N点的离散傅里叶变换 如果X的实际长度N,则X中的多余部分将被删除。快速傅立叶变换（FFT）是离散傅里叶变换的一种快速算法。（2）ifft函数：一维快速傅立叶逆变换 调用形式：Y=ifft（X，N）注：（1）当N取2的整数幂时，傅立叶变换的计算速度最快；通常取大于且最靠近X实际长度的幂次。N=2
33、p p=nextpow2（X的实际长度）（2）一般情况下，fft的结果为复数，可用abs及angle 分别求其幅度、相位。现在学习的是第32页，共40页例1、已知模拟信号求N=128点DFT的幅度谱和相位谱。MATLAB命令窗口输入如下：N=128;n=0:127;t=0:0.01:127;%时域数据点数 q=n*2*pi/N;%频域数据点数 X=2*sin(6*pi*t)+7*cos(8*pi*t);Y=fft(X,N);subplot(2,1,1);plot(q,abs(Y)subplot(2,1,2);plot(q,angle(Y)现在学习的是第33页，共40页例2、比较原序列和经过FF
34、T和IFFT变换后序列。MATLAB命令窗口输入如下：n=0:127;t=0:0.01:pi;%时域数据点数 q=n*2*pi/N;%频域数据点数 X=cos(80*pi*t);Y=real（ifft（fft(X)）;subplot(2,1,2);plot(t，X)title（原信号）subplot(2,1,2);plot(t，Y)title（两次变换后恢复的信号）现在学习的是第34页，共40页例3、利用傅立叶变换和卷积公式求两个离散序列的卷积。（1）利用卷积公式求卷积：f1=0，ones（1，9）；f2=0，0，0，ones（1，6）；f=conv（f1，f2）（2）利用傅立叶变换求卷积：N
35、=32 ；%傅立叶变换返回点数 F1=fft（f1，N）；F2=fft（f2，N）；F=F1.*F2；f=real（ifft（F）；现在学习的是第35页，共40页3.6数值积分数值积分Matlab提供了提供了quad、quadl、trapz等数值积分指令。等数值积分指令。quad：自适应：自适应Simpleson积分法计算数值定积分。积分法计算数值定积分。q=quad(fun,a,b)近似地从近似地从a到到b计算函数计算函数fun的数值积分，缺省误差为的数值积分，缺省误差为10-6。若给。若给fun输入向量输入向量x，应返回向量，应返回向量y，即，即fun是一单值函数。是一单值函数。q=qua
36、d(fun,a,b,tol)用指定的绝对误差用指定的绝对误差tol代替缺省误差。代替缺省误差。tol越大，函数计算的次数越少，速度越快，但结果精度变小。越大，函数计算的次数越少，速度越快，但结果精度变小。36现在学习的是第36页，共40页quadl:自适应自适应Lobatto积分法计算数值定积分，精度比积分法计算数值定积分，精度比quad高，指令格式与高，指令格式与quad相同。相同。fun=inline(3*x.2./(x.3-2*x.2+3);Q1=quad(fun,0,2)Q2=quadl(fun,0,2)trapz:梯形法计算数值积分梯形法计算数值积分，精度较低，但速度快。，精度较低，
37、但速度快。T=trapz(X,Y)用梯形法计算用梯形法计算Y在在X点上的积分。要求点上的积分。要求X和和Y为长度相同的向量。为长度相同的向量。X=-1:.1:1;Y=1./(1+25*X.2);T=trapz(X,Y)37现在学习的是第37页，共40页3.7常微分方程（常微分方程（ODE）组初值问题的数值解）组初值问题的数值解Matlab提供了提供了ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb等指令来解决常微分方程的数值解问题。这些指令在解等指令来解决常微分方程的数值解问题。这些指令在解题类型和精度上各有不同，其中在大多数场合的题类型和精度上各
38、有不同，其中在大多数场合的首选算法是首选算法是ode45.T,Y=ode45(odefun,tspan,y0)在区间在区间tspan=t0,tf上，从上，从t0到到tf，用初始条件，用初始条件y0求解显式微分方求解显式微分方程程y=f(t,y)。对于标量。对于标量t与列向量与列向量y，函数，函数f=odefun(t,y)必须返回一必须返回一f(t,y)的列向量的列向量f。解矩阵。解矩阵Y中的每一行对应于返回的时间列向量中的每一行对应于返回的时间列向量T中的一个时间点。要获得问题在其他指定时间点中的一个时间点。要获得问题在其他指定时间点t0,t1,t2,上的解，上的解，则令则令tspan=t0,
39、t1,t2,tf（要求是单调的）。（要求是单调的）。38现在学习的是第38页，共40页odefun为显式常微分方程为显式常微分方程y=f(t,y)，或为，或为包含一混合矩阵的方程包含一混合矩阵的方程M(t,y)*y=f(t,y)。命令命令ode23只能求解常数混合矩阵的问题；只能求解常数混合矩阵的问题；命令命令ode23t与与ode15s可以求解奇异矩阵的可以求解奇异矩阵的问题。问题。tspan:积分区间（即求解区间）的向量积分区间（即求解区间）的向量tspan=t0,tf。要获得问题在其他指定时。要获得问题在其他指定时间点间点t0,t1,t2,上的解，则令上的解，则令tspan=t0,t1,
40、t2,tf（要求是单调的）。（要求是单调的）。y0:包含初始条件的向量。包含初始条件的向量。39现在学习的是第39页，共40页例例求解描述振荡器的经典的求解描述振荡器的经典的VerderPol微分方程微分方程解：首先将该二阶方程化为一阶方程组。令解：首先将该二阶方程化为一阶方程组。令x1=y，x2=dy/dx，则，则dx1/dt=x2,dx2/dt=(1-x2)-x1编写函数文件编写函数文件verderpol.m：functionxprime=verderpol(t,x)globalMUxprime=x(2);MU*(1-x(1)2)*x(2)-x(1);再在命令窗口中执行：再在命令窗口中执行：globalMUMU=7;Y0=1;0t,x=ode45(verderpol,0,40,Y0);x1=x(:,1);x2=x(:,2);plot(t,x1,t,x2)40现在学习的是第40页，共40页