《数字信号处理》教学参考资料（MATLAB 5手册）第11章 积分和微分方程组

China-pub.com 下载 第11章 积分和微分方程组 在有效的MATLAB命令帮助下，可以求解出定积分和普通微分方程的数字解并绘制出其 图形。 11.1 积分 在MATLAB中能求解如下形式的定积分并给出数字解： g=g(x)dx 有许多方法都可以能够解决积分问题（又叫做求面积）。如果要用MATLAB监控整个计算 过程，可以使用guad命令。同样能计算出被积函数g的值，并且让MATLAB使用梯形规则和 trapz命令计算出积分。当只有离散的数据点和被积函数的数学表达式为未知时，这种方法 是非常有效的。 命令集107 定积分计算 trapz(x,y) 计算出函数的积分并将结果返回到。向量x和y有相同的长度， i,y代表曲线上的一点。曲线上点的距离不一定相等，x值也 不一定有序。然而，负值间距和子区间被认为是负值积分。 trapz(y) 计算方法同上，但x值间隔为1。 trapz(x,A) 将A中每列的值带入x的函数算出其积分，并返回一组包含 积分结果的向量。A的列向量必须和向量x的长度相同。 Z=trapz(x,A,dim) 在矩阵A中dim指定的维内进行数据积分。如果给定向量x, 则x的长度必须与size(A,dim)相同。 cumtrapz(A,dim) 返回大小和A相同的数组，包含的是将矩阵A进行梯形积分 的累积值。如果dim已给定，则在dim维内进行计算。 quad(fcn,a,b) 返回在区间a,b]上g的积分近似值。字符串fcn包含一个与g相 对应的MATLAB函数名，也就是预定义函数或者是M文件。 这个函数接收一个向量参数，并返回一个向量结果。 MATLAB利用辛普森规则执行递归的积分，计算误差为1O~。 guad(fcn,a,b,tol)求g的积分近似值，其相对误差由参数tol定义。否则，计算 过程同上。 quad(fcn,ab,tol求g的积分近似值，其相对误差由参数tol所定义。如果参数 pic) pic是非零值，则在图形中显示求值的点。 quad(···,trace) 如果trace是非零值，则画出积分图形

下载 第11章 积分和微分方程组 在有效的M AT L A B命令帮助下，可以求解出定积分和普通微分方程的数字解并绘制出其 图形。 11.1 积分 在M AT L A B中能求解如下形式的定积分并给出数字解： 有许多方法都可以能够解决积分问题 (又叫做求面积 )。如果要用 M AT L A B监控整个计算 过程，可以使用 q u a d命令。同样能计算出被积函数 g的值，并且让M AT L A B使用梯形规则和 t r a p z命令计算出积分。当只有离散的数据点和被积函数的数学表达式为未知时，这种方法 是非常有效的。 命令集1 0 7 定积分计算 t r a p z ( x , y ) 计算出函数x的积分并将结果返回到y。向量x和y有相同的长度， (xi, yi)代表曲线上的一点。曲线上点的距离不一定相等，x值也 不一定有序。然而，负值间距和子区间被认为是负值积分。 t r a p z ( y ) 计算方法同上，但x值间隔为1。 t r a p z ( x , A ) 将A中每列的值带入 x的函数算出其积分，并返回一组包含 积分结果的向量。A的列向量必须和向量x的长度相同。 Z = t r a p z ( x , A , d i m ) 在矩阵A中d i m指定的维内进行数据积分。如果给定向量 x， 则x的长度必须与size(A, dim)相同。 c u m t r a p z ( A , d i m ) 返回大小和A相同的数组，包含的是将矩阵 A进行梯形积分 的累积值。如果d i m已给定，则在d i m维内进行计算。 q u a d ( f c n , a , b ) 返回在区间[a, b]上g的积分近似值。字符串f c n包含一个与g相 对应的M AT L A B函数名，也就是预定义函数或者是 M文件。 这个函数接收一个向量参数，并返回一个向量结果。 M AT L A B利用辛普森规则执行递归的积分，计算误差为1 0－3。 q u a d ( f c n , a , b , t o l ) 求g的积分近似值，其相对误差由参数 t o l定义。否则，计算 过程同上。 quad(fcn,a b,tol,求g的积分近似值，其相对误差由参数 t o l所定义。如果参数 p i c ) p i c是非零值，则在图形中显示求值的点。 q u a d (. . ., t r a c e ) 如果t r a c e是非零值，则画出积分图形

152 China-pub.com MATLAB5手册 下载 guad8(···) 可以与quad一样用于相同的参数组合并返回相同的结果， 但使用更高精度的方法。因此，如果被积函数的导数在某一 区间内是不定的，例如：9=5md,使用此命令将会更好 一些。quad和quad8都要求被积函数在整个区间里是有限的。 dblquad(f,min1,计算双变量函数f的二重积分。函数中的第一个自变量用于 max1,min2,max2, 内层积分。内层积分在minl和nax 1之间进行，外层积分在 tol,trace,order) min2和max2之间进行。变量tol指定相对误差。trace的使用 方法与quad相同。根据字符串order,对于相同的访问， dblquad能选择使用quad、quad8和许多用户定义的积分 方法，并返回与guad相同的变量。 输入quaddemo可以看到一个演示实例。 ■例11.1 下面用不同的方法来计算下列积分： (a)使用trapzt命令。首先创建一个有x值的向量。用5和l0两个值进行计算： x5=1 inspace(0,1,5);x10=1 inspace(0,1,10); 然后创建x的函数y: y5=exp(-x5.2);y10=exp(-x10.2); 现在计算出积分值： format long; integral5 trapz(x5,y5),... integral10 trapz(x10,y10) 返回 integral5 0.74298409780038 integral10 0.74606686791267 (b)使用guad命令。首先在M文件中创建函数。此文件integrand.m包含函数，如下： function y integrand(x) y =exp(-x.-2); 首先以标准误差计算积分，然后再以指定误差计算积分。 format long; integralStd quad('integrand',0,1) integralTol quad('integrand',0,1,0.00001)

q u a d 8 (. . .) 可以与q u a d一样用于相同的参数组合并返回相同的结果， 但使用更高精度的方法。因此，如果被积函数的导数在某一 区间内是不定的，例如： ,使用此命令将会更好 一些。q u a d和quad8 都要求被积函数在整个区间里是有限的。 dblquad(f, min1,计算双变量函数f的二重积分。函数中的第一个自变量用于 m a x 1 , m i n 2 , m a x 2 , 内层积分。内层积分在m i n1 和m a x1之间进行，外层积分在 t o l , t r a c e , o r d e r ) m i n2和m a x2之间进行。变量t o l指定相对误差。 t r a c e的使用 方法与 q u a d相同。根据字符串 o r d e r，对于相同的访问， d b l q u a d能选择使用q u a d、q u a d 8和许多用户定义的积分 方法，并返回与q u a d相同的变量。 输入q u a d d e m o可以看到一个演示实例。 ■ 例11 . 1 下面用不同的方法来计算下列积分： (a) 使用t r a p z命令。首先创建一个有x值的向量。用5和1 0两个值进行计算： 然后创建x的函数y: 现在计算出积分值： (b) 使用q u a d命令。首先在M文件中创建函数。此文件i n t e g r a n d . m包含函数，如下： 首先以标准误差计算积分，然后再以指定误差计算积分。 q = sinx 0 1 ò dx 1 5 2 M ATLAB 5 手册 下载 返回

China-pub.com 第11章积分和微分方程组 153 下载 给出 integralStd 0.74682612052747 integralTol 0.74682414517798 (c)使用quad8命令：使用在(b)中创建的M文件，然后输入： integral8Std quad8('integrand',0,1) integral8Std 0.74682413281243 这是MATLAB所能给出的最精确的结果。 (d使用cumtrapz命令能很容易地计算出不同区间的积分。 x=0:5; cumtrapz(x) ans 00.5000 2.0000 4.5000 8.0000 12.5000 (e)计算二重积分： 如图11-l。首先创建一个包含函数M文件：integrand.2.m: function y integrand2(x,y) y=exp(-x.2-y.2); 然后用quad命令计算对于固定的x值在y方向的一些积分值： x=linspace(0,1,15); for1=1:15 integral(i)=quad('integrand2',0,1,[],[],x(i)) end 现在已计算出在y方向的15个积分值。trapz命令能使用这些值来计算二重积分： format short; dIntegral trapz(x,integral) dIntegral 0.5575 输入下列语句可以得到一个积分区域的图形： [X,Y]=meshgrid(0:.1:1,0:.1:1); Z =integrand2(X,Y); mesh(X,Y,Z);view(30,30);

(c) 使用q u a d 8命令：使用在( b )中创建的M文件,然后输入： 这是M AT L A B所能给出的最精确的结果。 (d) 使用c u m t r a p z命令能很容易地计算出不同区间的积分。 (e) 计算二重积分： 如图11 - 1。首先创建一个包含函数M文件:i n t e g r a n d 2 . m: 然后用q u a d命令计算对于固定的x值在y方向的一些积分值： 现在已计算出在y方向的1 5个积分值。t r a p z命令能使用这些值来计算二重积分： 输入下列语句可以得到一个积分区域的图形： 第11章 积分和微分方程组 1 5 3 下载 给出

154 China-pub.coM MATLAB5手册 下载 结果如图11-1所示。命令mesh和view定义在13.5节中。 0.日 0. 。0 图11-1函数e2在区间[0,1]×[0,1]的上的图形 ■ 不定积分 ∫f)d不能使用上面的命令来计算。MATLAB中的数学符号工具箱和 MATLAB的编辑器能提供处理这些积分的命令。 11.2常微分方程组 下面来研究常微分方程系统ODE,该系统处理的是初始值已知的一阶微分方程。在本节中 主要讨论这种类型的微分方程，同时也会举出两个有关边界值问题的例子。可以利用ODE系 统创建稀疏线性系统方程来求解这些例子。 在数学符号工具箱中，有一些命令能给出常微分方程的符号解，即解以数学表达式的形式给出。 在下面的初始值问题中，有两个未知函数：x1()和x2(),并用以下式子表达其微分形式： d水i二x对 d 在许多应用中，独立变量参数表示时间。 x=f(x1,x2,t) x2=f2(x1,x2,t) x1(0)=x1,0 x2(0)=x2,0 高阶的ODE能表达成第1阶的ODE系统。例如，有以下微分方程： x”=f(x,x',t) x(0)=x0 x'(to)=xpo 用x替换x‘用x,替换x,就能得到： [x1=2 x2=f(x1,x2,t) x1(0)=x0 x2(0)=xp0

结果如图11 - 1所示。命令m e s h和v i e w定义在1 3 . 5节中。 图11-1 函数e -x2 -y2在区间[ 0，1 ]×[ 0，1 ]的上的图形 不定积分 不能使用上面的命令来计算。 M AT L A B中的数学符号工具箱和 M AT L A B的编辑器能提供处理这些积分的命令。 11.2 常微分方程组 下面来研究常微分方程系统 O D E, 该系统处理的是初始值已知的一阶微分方程。在本节中 主要讨论这种类型的微分方程，同时也会举出两个有关边界值问题的例子。可以利用 O D E系 统创建稀疏线性系统方程来求解这些例子。 在数学符号工具箱中，有一些命令能给出常微分方程的符号解，即解以数学表达式的形式给出。 在下面的初始值问题中，有两个未知函数：x 1(t)和x 2(t)，并用以下式子表达其微分形式： 在许多应用中，独立变量参数 t表示时间。 高阶的O D E能表达成第1阶的O D E系统。例如，有以下微分方程： 用x2替换x 用x1替换 x，就能得到： f(t) a x ò dt 1 5 4 M ATLAB 5 手册 下载 dxi dt = x ¢ i ■

China-pub.coM 第11章积分和微分方程组 155 下载 这是一个第1阶的ODE系统。 对于某一时间间隔0≤1≤T,初始值问题的解决方法是将时间分成一组有限和离散的时间 点，例如用相同的时间间隔△进行等分： t=i△t,i=0,.,W 其中时间步长△仁TN,N为某一整数。这种导数能被微分方程的可微分的商所代替，微分 方程表示在不同时间点的解。见例11.2，给出更多的有关有限微分商的信息。这种方法的稳 定性取决于△的大小和所采用的数值方法，用这种方法能得到ODE的近似值。 在许多应用中有一些微分过程非常复杂的微分方程，在某些区域里这些方程要求有非常 小的时间步长△1。解决这些问题的困难在于问题中涉及不同的时间尺度，如解的导数可能有 较大的变化。 MATLAB使用龙格-库塔-芬尔格(Runge-Kutta-Fehlberg)方法来解ODE问题。在有限点内 计算求解，而这些点的间距由解本身来决定。当解比较平滑时，区间内使用的点数少一些： 在解变化很快时，区间内应使用较多的点。 为了得到更多的有关何时使用哪种解法和算法的信息，推荐使用helpdesk。所有求解 方程通用的语法或句法在命令集108中头两行给出。时间间隔将以向量t=[t0,tt]给出。 命令ode23可以求解(2,3)阶的常微分方程组，函数ode45使用(4,5)阶的龙格-库塔-芬 尔格方法。注意，在这种情况下x是x的微分，不是x的转置。 在命令集108中so1ver将被诸如ode45函数所代替。 命令集108 龙格-库塔-芬尔格方法 [time,x]= 计算ODE或由字符串str给定的ODE的值。部分解已 solver(str,t,x0) 在向量time中给出。在向量time中给出部分解，包含的 是时间值。还有部分解在矩阵X中给出，X的列向量是 每个方程在这些值下的解。对于标量问题，方程的解将 在向量X中给出。这些解在时间区间(1)到t(2)上计算得 到，其初始值是x0即x(t(1)。此方程组由str指定的M文 件中函数表示出。这个函数需要两个参数：标量1和向 量x,应该返回向量x'(即x的导数)。因为对标量ODE来 说，x和x'都是标量。在M文件中输入odefile可得到 更多信息。同时可以用命令numjac来计算雅可比函 数。 [t,X]= 此方程的求解过程同上。结构val包含用户给solver solver(str,t,x0,val) 的命令。参见odeset和表11-1，可得更多信息。 ode45 此方法被推荐为首选方法。 ode23 这是一个比ode45低阶的方法。 ode113 用于更高阶或大的标量计算。 ode23t 用于解决难度适中的问题。 ode23s 用于解决难度较大的微分方程组。对于系统中存在常 量矩阵的情况也有用

这是一个第1阶的O D E系统。 对于某一时间间隔0≤t≤T，初始值问题的解决方法是将时间分成一组有限和离散的时间 点，例如用相同的时间间隔 Dt进行等分： 其中时间步长Dt=T / N，N为某一整数。这种导数能被微分方程的可微分的商所代替，微分 方程表示在不同时间点的解。见例 11 . 2，给出更多的有关有限微分商的信息。这种方法的稳 定性取决于Dt的大小和所采用的数值方法，用这种方法能得到 O D E的近似值。 在许多应用中有一些微分过程非常复杂的微分方程，在某些区域里这些方程要求有非常 小的时间步长 Dt。解决这些问题的困难在于问题中涉及不同的时间尺度，如解的导数可能有 较大的变化。 M AT L A B使用龙格-库塔-芬尔格(R u n g e - K u t t a - F e h l b e rg)方法来解O D E问题。在有限点内 计算求解，而这些点的间距由解本身来决定。当解比较平滑时，区间内使用的点数少一些； 在解变化很快时，区间内应使用较多的点。 为了得到更多的有关何时使用哪种解法和算法的信息，推荐使用 h e l p d e s k。所有求解 方程通用的语法或句法在命令集 1 0 8中头两行给出。时间间隔将以向量 t = [ t 0 , t t ]给出。 命令o d e 2 3可以求解( 2，3 )阶的常微分方程组，函数 o d e 4 5使用( 4，5 )阶的龙格-库塔-芬 尔格方法。注意，在这种情况下 x´是x的微分，不是x的转置。 在命令集1 0 8中s o l v e r将被诸如o d e 4 5函数所代替。 命令集1 0 8 龙格-库塔-芬尔格方法 [ t i m e , X ]= 计算O D E或由字符串s t r给定的O D E的值。部分解已 s o l v e r ( s t r , t , x 0 ) 在向量t i m e中给出。在向量t i m e中给出部分解，包含的 是时间值。还有部分解在矩阵 X中给出，X的列向量是 每个方程在这些值下的解。对于标量问题，方程的解将 在向量X中给出。这些解在时间区间t ( 1 )到t ( 2 )上计算得 到，其初始值是x 0即x ( t ( 1 ) )。此方程组由s t r指定的M文 件中函数表示出。这个函数需要两个参数：标量 t和向 量x，应该返回向量x' (即x的导数)。因为对标量O D E来 说，x和x'都是标量。在M文件中输入o d e f i l e可得到 更多信息。同时可以用命令 n u m j a c来计算雅可比函 数。 [ t , X ] = 此方程的求解过程同上。结构 v a l包含用户给s o l v e r s o l v e r ( s t r , t , x 0 , v a l ) 的命令。参见o d e s e t和表11 - 1，可得更多信息。 o d e 4 5 此方法被推荐为首选方法。 o d e 2 3 这是一个比o d e 4 5低阶的方法。 o d e 1 1 3 用于更高阶或大的标量计算。 o d e 2 3 t 用于解决难度适中的问题。 o d e 2 3 s 用于解决难度较大的微分方程组。对于系统中存在常 量矩阵的情况也有用。 第11章 积分和微分方程组 1 5 5 下载