《精通MATLAB科学计算（第3版）(王正林)12-3r.pdf》由会员分享，可在线阅读，更多相关《精通MATLAB科学计算（第3版）(王正林)12-3r.pdf（32页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第 章常微分方程求解常微分方程作为微分方程的基本类型之一，在自然界与工程界有很广泛的应用。很多问题的数学表述都可以归结为常微分方程的定解问题。很多偏微分方程问题，也可以化为常微分方程问题来近似求解。MATLAB中提供了求解函数和求解器，可以实现常微分方程的解析求解和数值求解，而且还可以通过编程实现一些常见的数值解法。通过本章学习，读者能够熟练使用MATLAB的求解函数和求解器，以及通过编程进行常微分方程的求解。12.1MATLAB中的求解函数一般微分方程式描述系统内部变量的变化率如何受系统内部变量和外部激励，如输入的影响。当常微分方程式能够解析求解时，可用MATLAB的符号工具箱中的功能找到精
2、确解;在常微分方程难以获得解析解的情况下使用MATLAB的常微分方程求解器solver,可以方便地在数值上求解。12.1.1 符号解函数dsolve一阶常微分方程式(first-order Ordinary Differential Equation,ODE)可写为：y=g(x,y)其中x 为独立变量，而)，是 x 的函数。它的解是：y=.f(x,y)该解可以满足)/=f=g(x j),在初始条件.i，(xo)=yo下，可以得到唯一解。MATLAB常微分方程符号解的语法是：dsolve(equation*,1 condition 1)其 中，equation代表常微分方程式即./=g Q y)
3、,且须以Dy代表一阶微分项y ,D2y代表二阶微分项y,condition则为初始条件。函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。第12章常微分方程求解dsolve的调用格式如下：(1)dsolve(equation)给出微分方程的解析解，表示为，的函数；(2)dsolve(equation condition)给出微分方程初值问题的解，表示为/的函数；(3)dsolveCequation,vf)给出微分方程的解析解，表示为v 的函数；(4 )dsolve(*equation,condition*,V)给出微分方程初值问题的解，表示
4、为v 的函数。【例 12-1】常微分方程符号解求解实例1。计算微分方程或+3孙=配-的通解。dx dsolve(Dy+3*x*y=x*exp(-xA2)1)输出为：ans=(l/3*exp(-x*(x-3*-)+C1)*exp(-3*x*t)t)+C1)*士 3*、&)由于系统默认的自变量是/,显然系统把x 当作常数，把y 当作/的函数求解.输入命令:dsolve(Dy+3*x*y=x*exp(-xA2)z*x)输出为：ans=exp(-x、)+exp(-3/2*xA2)*C1exp(-xA2)+oxp(3/2*M 2)*C1_3 2可知通解为y=*G，其中G 为常数。【例 12-2】常微分方
5、程符号解求解实例2。计算微分方程中+2 y-/=0 在初始条件川”1=2e下的特解。dsolve(*x*Dy+2*y-exp(x)=0,*y(1)=2*exp(1),*x)输出为：ans=(xp(x)*x-exp(x)+2*xxp(1)/x八2(oxp(x)*x-o xp (x)+2*派断(44-)-笠可知特解为y=+2 e。X【例 12-3】常微分方程符号解求解实例3。求炉+2丁+产=0 的通解。dsolve(*D2y+2*Dy+exp(x)=0,1x*)Y 281精通MATLAB科学计算(第2忡-输出为：ans=-l/3*exp(x)-l/2*exp(-2*x)*C1+C2可知通解为,-L
6、 N*。？，其中G,C2为常数。3 2282 第12章常微分方程求解12.1.2 求解器 solver在求常微分方程数值解方面，MATLAB具有丰富的函数，将其统称为solver,其一般格式为：Tf Y=solver(odefun,tspan,yO)该函数表示在区间tspan=Do,5上,用初始条件川求解显式常微分方程V=odefun为显式常微分方程V=f(t,y)中的/(f,y),tspan为求解区间，要获得问题在其他指定点2上的解,则令他劭=汨/1,/2,0(要 求 单 调)，y 0 为初始条件。solver 为命令 ode45、ode23、odel 13,ode 15s,ode23s、o
7、de23t、ode23tb 之 一，其中ode45、ode23、odel 13属于非刚性ODE类 型,这些命令的特点如表12-1所 示;ode 15s,ode23s、ode23t、ode23tb属于刚性ODE类 型，这些命令的特点如表12-2所示。表 1 2-1 非刚性O D E 求解命令求解器s o l v e r功 能说 明o d e 45一步算法；4,5 阶龙格库塔方程；累计截断误差达(&03大部分场合的首选算法o d e 23一步算法；2,3 阶龙格-库塔方程；累计截断误差达(加个使用于精度较低的情形o d e l 13多步法；Ad a m s 算 法；高低精度均可到IO。I。/,计算
8、时间比o d e 45短表 1 2-2 刚性O D E 求解命令求解器s o l v e r功 能说 明o d e 23t梯形算法适度刚性情形o d e 15s多步法；G e a r s 反向数值微分；精度中等若 o d e 45失效时，可尝试使用o d e 23s一步法；2 阶 Ro s e b r o c k 算 法；低精度当精度较低时，计算时间比o d e l 5s 短o d e 23t b梯形算法；低精度当精度较低时，计算时间比o d e l 5s 短函数ode45与 ode23是常使用的求解方法，函数ode45的使用与ode23完全一样。两个函数的差别在于必须与所用的内部算法相关。两
9、个函数都运用了基本的龙格-库塔(Runge-Kutta)数值积分法的变形。ode23运用一 组 合 的 2/3阶龙格-库塔-芬尔格(Runge-Kutta-Fehlerg)算 法，而 ode45运用组合的4/5阶龙格-库塔-芬尔格算法。通 常，ode45可取较多的时间步。所 以，要保持与ode23相同误差时，在访和夕之间可取较少的时间步。然 而，在同一时间内，ode23每时间步至少调用3 次，而 ode45每时间步至少调用6 次。正如使用高阶多项式内插常常得不到最好的结果一样，ode45也不总是比。de23好。如果。de45产生的结果，对作图间隔太大，则必须在更细的时间区间内对数据进行内插，比
10、如用函数interpl。这个附加时间点会使ode23更有效。作为一条普遍规则，在所计算的1 283精通MATLAB科学计算（第2版i-导数中，如有重复的不连续点，为保持精度致使高阶算法减少时间步长，这时低阶算法更有效。采用求解器solver求解ODE的基本过程如下：（1 ）根据问题所属学科中的规律、定律、公 式，用微分方程与初始条件进行描述。F（y,y,.yn-t）=0y（0）=Jo,y（O）=i,.yn-l（0）=y-写为向量的形式为y=y；MD；M2）,丁（m一1）,与加可以不等。（2）运用数学中的变量替换：为=y（n-l）,为/可（n2）,.,yi=y=y，把高阶（大于2阶）的方程（组）
11、写成一阶微分方程组：%1y=力”,y)力(“)7相应的初始条件为：（3）根 据（1 ）与（2）的结果，编写能计算导数的M-函数文件odefileo（4）将文件odefile与初始条件传递给求解器solver中的一个，运行后就可得到ODE在指定时间区间上的解列向量y（其中包含j，及不同阶的导数【例 12-4 求解器solver应用实例。采用ode45求解描述某非刚性体的运动方程的”=y2y3 伊（o）=o微分方程：yi=一川3，其初始条件为了 2（。）=1 ,并画出解的图形。y 3=0.5112 乃（。）=1解：编写函数文件rigid.m:function dy=rigid(tz y)dy=ze
12、ros(3,1);%行向量dy(1)=y(2)*y(3);dy(2)=-y(l)*y(3);dy(3)=-0.51*y(l)*y(2);建立文件exl204.m,如下所示：284 第12章常微分方程求解%exl204.moptions=odeset(*RelTol*zle-4,1AbsTol1,le-4 le-4 le-5);T,Y=ode45(rigid,0 12,0 1 1,options);plot(TZY(:Z1)，,T,Y(:,2)J-1T,Y(:,3),1 J)legend(yl y2 y3)运行程序，输出的图形结果为图12-1所示。12.2欧拉法图 1 2-1 非刚性体运动的微分
13、方程图12.2.1简单欧拉法欧拉法(Euler)是简单有效的常微分方程数值解法，欧拉法有多种形式的算法，其中简单欧拉法是一种单步递推算法。1.基本原理对于一个简单的一阶微分方程：y=f(x,y)其中初始值为：y(x0)=y0欧拉方法等同于将函数微分转换为数值差分，如下式,y(x)J.+h故原方程变形为：yn+=y+hf(x,l,y)i 285精通MATLAB科学计算(第2版i-2.算法的程序实现在M ATLAB中编程实现的欧拉法函数为：MyEulero功 能：以欧拉法求解常微分方程。调用格式:outx.outy=MyEuler(fun,xO,xt,jX),PointNum)0其 中，fun为函
14、数名；xO为自变量区间初值；“为自变量区间终值；可为函数在xO处 的 值；PointNum为自变量在M),xf上所取的点数；outx为输出自变量区间上所取点的x值；outy为对应点上的函数值。欧拉法的M ATLAB程序代码如下：f u n c t i o n o u t x,o u t y=My Eu l e r(f u n,x Oz x t,y O,P o i n t Nu m)金用前向差分的欧拉方法解微分方程y =f u n率函数 f (x,y)：f u n务自变量的初值和终值：x Ozx t%y 0表示函数在x O处的 值，输入初值为列向量形式 自变量在 x O,x t 上取的点数:P
15、o i n t Nu m%o u t x:所取的点的x值%o u t y:对应点上的函数值i f n a r g i n 5|P o i n t Nu m=0P o i n t Nu m=100;e n di f n a r g i n+i=yn+5/,y”)+/(x+i,yn+)2.算法的程序实现在 MATLAB中编程实现的改进的欧拉法函数为：MyEulerProo功 能：用改进的欧拉法求解常微分方程。调用格式:Xout,Yout=MyEulerPro(fun,xO,xt,PointNumber)o其 中，ftin为函数名；M)为自变量区间初值；X，为自变量区间终值；川为函数在xO处 的 值
16、；288 第 1 2 章常微分方程求解P o i n t Nu m b e r 为自变量在 x O,词上所取的点数；Xo u t 为输出自变量区间上所取点的x值；Yo u t 为对应点上的函数值。改进的欧拉法的M A T L A B 程序代码如下：function Xout,Yout=MyEulerPro(fun,xO,xtz yO,PointNumber)%用改进的欧拉法解微分方程y，=fun%函数 f(x,y)：fun当自变量的初值和终值：xO,xt%y0表示函数在xO处的 值，输入初值为列向量形式许自变量在 xO,xt 上取的点数:PointNumber%Xout:所取的点的x 值%Yo
17、ut:对应点上的函数值if nargin5 I PointNumer=0为 100PointNumer=100;endif nargin4%y0 默认值为 0y0=0;endh=(xt-xO)/PointNumber;x=x0+0:PointNumber *h;y(lr:)=y0(:);for i=l:PointNumberfl=h*feval(fun,x(i),y(i,:告如果函数仅输入4 个参数值，则 PointNumer默认值为计算所取的两离散点之间的距离%表示出离散的自变量x常迭代计算过程)；fl=fl(:)f2=h*feval(fun,x(i+1)zy(i,:)+fl);f2=f2(
18、:)y(i+l,:)=y(i,:)+1/2*(fl+f2);endXout=x;Yout=y;3.实例分析【例 12-6】改进的欧拉法求解常微分方程应用实例。利用改进的欧拉法求常微分方程：y=s i n x +yy(x o)=l,x o =0即 对【例 12-4用改进的欧拉法求解，比较两种解法的结果，并画出解的图形。解：编写 M A T L A B 程序 e x l 206.m ：function exl206()289精通MATLAB科 学 计 算(第2蝌图 1 2-3 改进的欧拉法解、欧拉解、精确解的对比图和 e x l 206比较改进的欧拉法,简单欧拉方法以及微分方程符号解 x 3,y
19、3=My Eu l e r P r o(*m y f u n 01,0,2*p i,1,128);x,y l =My Eu l e r (*m y f u n 01*,0,2*p i,1,128);超欧拉法所得的解y=d s o l v e (Dy=y+s i n (t),y (0)=11);超该常微分方程的符号解f o r k=l:129t(k)=x(k);y 2(k)=s u b s(y,t(k);招求其对应点的离散解e n dp l o t (x,y l J -b l x 3,y 3,Pg,x,y 2 J *r，)l e g e n d (，简单欧拉法解，J 改进欧拉解，J 符号解，)运
20、行程序后，输出如图12-3所示的图形。如 图 12-3所 示，改进欧拉解与精确解几乎完全吻合，而简单欧拉解与精确解还有一定的误差。由 此 可 见，改进的欧拉法较之简单欧拉法更为精确。290 第12章常微分方程求解12.3龙格甚塔法尽管改进的欧拉法相对于简单欧拉法较为精确。但是对于很多实际的问题，运用这两种方法仍然得不到足够精确的解。龙格-库塔法(Runge-Kutta)是较之前两种方法计算精度更高的方法。1.基本原理在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为。()，被广泛应用于解微分方程的初值问题。其算法公式为：h%+1 yn+二(ki+2k2+2k3+左4 )6其 中：ki=f(x
21、,yn)ki=Jxn+h,y+；/?肩)左 3 =f(x+；万 左 2)k4=f(x+h,y+府3)2.四阶龙格-库塔算法编程实现在 MATLAB中编程实现的四阶龙格-库塔算法函数为：MyRunge_Kuttao功 能：用四阶龙格-库塔算法求解常微分方程。调用格式:x,y=MyRunge_Kutta(ftin,xO,xl,yO,PointNum,varargin)o其 中，fbn为函数名；xO为自变量区间初值；口为自变量区间终值；yO为函数在M)处的值；PointNum为自变量在xO,xf上所取的点数；Varargin为可选项参数；x 为输出自变量区间上所取点的x 值；y 为对应点上的函数值。
22、四阶龙格-库塔法的MATLAB程序代码如下：fu n c tio n x,y=M yRunge_Kutta(fun,xOz x tz yO,P ointN um,varargin)%R unge-K utta方法解微分方程形为y(t)=f(x,y(x)考此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式 291精通MATLAB科 学 计 算(第2蚪殳函数 f(x,y)：fun,自变量的初值和终值：xO,xtSyO表示函数在xO处的 值，输入初值为列向量形式*自变量在 xO,xt上取的点数:PointNum%varargin为可选输入项可传适当参数给函数f(x,y)%x:所取的点的
23、x值务y:对应点上的函数值if nargin 4|PointNum=0PointNum=100;endif nargin 第12章常微分方程求解常用三种不同方法解微分方程clear,elfxO=O;xt=2;Num=100;h=(xt-xO)/(Num-1);x=xO+0:Num*h;a=1;yt=1-exp(-a*x);fun=inline(-y+l*z*x*z*y);yO=0;PointNum=4;告清内存中的变量告真值解招用inline构造函数f(x,y)号设定函数初值%设定取点数xl,ye =MyEuler(fun,xO,xt,y0f PointNum);x2z yh=MyEulerP
24、ro(fun,xO,xt,y0z PointNum);x3f yr =MyRunge_Kutta(fun,xO,xtz yO,PointNum);plot(xz ytz 1k z xlz ye,b:1,x2,yh,*g:,x3z yr,r:)legend(1真值1,1简单欧拉法解1,改进的欧拉法解r,*Runge_Kutta法解1)hold onplot(xl,ye,*bo*,x2,yh,*b+,x3,yr,1r*1)PointNum=1000;%估计各算法运行PointNum步迭代的时间t i c电计时开始xl,ye =MyEuler(fun,xO,xt,yO,PointNum);time_
25、Euler=toe 带得到欧拉法的运行时间ticxlz yh=MyEulerPro(fun,xO,xtz yO,PointNum);time_EulerPro=toe 专 得到改进的欧拉法运行时间ticxlz yr=MyRunge_Kutta(fun,x0z xtz yO,PointNum);time_RK4=toe 宅 得至！I Runge_Kutta法运彳亍时间运行后得到：exl207time_Euler=0.25440.2544time_EulerPro=0.48870.4887time_RK4=1.01951.0195运行的结果如图12-4所示。Y 293精通MATLAB科 学 计 算
26、(第2蝌由运算结果可知，龙格-库塔法可得到与真值几乎相同的微分方程数值解。但其相应的程序运行时间较长，几乎是简单欧拉法的五倍(欧拉法运行时间是0.2544s,而龙格-库塔法运行时间是1.0195s卜4.求解器solvei中的龙格-库塔法求解前面讲到，求解器solver中的ode23采用二阶、三阶龙格-库塔法；ode45采用四阶、五阶龙格-库塔法。【例12-8】求解器solver中的龙格-库塔法求解应用实例1。分别采用二、三、四、五阶龙格-库塔方法求解以下方程，y=-3 y2+2x2+3x,H 0 x 1j(x o)=l,xo=0解法一：用二阶、三阶龙格-库塔函数求解。编写程序文件ex 1208
27、a.m:%exl208a.m用ode23得到微分方程解并计算出该算法运行时间fun=inline(-3*yA2+2*x.A2+3*x,x,y);%用 inline 构造函数 f(x,y)x,y=ode23(fun,(0,1,1);3可得到 x,y 输出向量值ode23(fun,0,1,1),hold on&可得到输出的函数图结果如图12-5a所不。294 -第 章常微分方程求解解法二：用四阶、五阶龙格-库塔函数求解。编写程序文件ex 1208b.m:%exl208b.m用ode45得到微分方程解，并计算出该算法运行时间fun=inline(-3*yA2+2*x.A2+3*x x1,*yT);告
28、 用 inline 构造函数 f(x,y)ode45(fun,0,1,1),hold on%可得到输出的函数图ticx,y=ode4 5(fun,0,1,1);tl=tocticx,y=ode23(fun,0,1,1);t2=toc命令窗口显不tl、t2 值。exl208btl=0.02560.0256t2=0.01660.0166如 图 12-5b所示显示函数的解结果。295精通MATLAB科 学 计 算(第2蝌由图12-5a和 图 12-5b所示可知，用两个函数求解得到的结果相同。而 ode23比 ode45的运行速度要快一些。【例 12-9】求解器solver中的龙格-库塔法求解应用实例
29、2。采用ode45求解如下二阶方程：V(2)=X2)-0.5+0.02XD初始条件为：x=0 时，XI)=25,y(2)=2,并画出解的图形。解：用 ode45求解。首先建立函数文件myfun02.m:myfun02.mfunction dx=myfun02(xz y)dx=zeros(2,1);dx(l)=y(l)*(l-0.1*y(2);dx(2)=y(2)*(-0.5+0.02*y(1);编写exl209.m对微分方程求解%exl209.mx,y=ode45(myfun02z 0 15,25 2);plot(x,y(:,l),-,x,y(:,2),*),画出 y(1),y(2)的函数图运
30、行程序，输出结果如图12-6所示。296 常微分方程求解12.4预估-校正法所谓预估-校正(p r e d i c t o r-c o r r e c t o r )算法是指算法中主要包括的两个过程，预估过程和校正过程。前面讲到的改进的欧拉法是一种简单的预估-校正算法。此外还有两种比较常用的形式，它们分别是 A B M (Ad a m s-Ba s h f b r t h-Mo u l t o n )法和 H a m m i n g 法。12.4.1 ABM 法A B M法是最常用的传统的预估-校正法。1.基本原理A B M法与上面几种方法的主要不同是：它属于多步算法，%的值是由，,(X._
31、2,y”-2),(x“_ 3,%-3)这几个点通过预估-校正过程求得的，要分以下两步进行。第一步:用四阶的拉格朗日多项式来近似得到/(XJ)的积分，并用前四个点的/值表示。从而可得y a+i的预估表达式p“+i ：pn+=yn+&(-9加3+37加2-59/t +55,/)第二步：对p向的预估值进行校正，得到为+i的校正值金+|。c“+i =%+一 5T+19人+9/m)通常使用的是改进的A B M算法。即在上述迭代式的基础上，再加入更正公式，因此可得如下的迭代运算公式：297精通MATLAB科学计算(第 2 忡-P+1=y.+&(-9 加 3+37小-59 小 +556)24加 +1 251
32、 z、Pn+Pn)。+1 =yn+(/一 2 一 5力】一 1 +19/+9/(X+,”+)19,、y/i+i=c+i 270 C，+1-P z)由于算法采用多步近似，因而其具有较小的截断误差。5),其计算精度优于前面介绍的众多方法。2.算法的程序实现MATLAB中的odel 13函数就是该算法的程序实现，由于迭代过程较复杂，此处不给出具体程序。有兴趣的读者可以参考odel 13的实现程序。其具体命令格式与上一节介绍的 ode23和 ode45相同。12.4.2 Hamming 法Hamming法是另一种形式的多步预估一校正法，其截断误差也是。(二),运算效率与ABM相当。下面对其进行具体介绍
33、。1.基本原理Hamming预估-校正公式为：4h“Pn+=y-3+(2.AI-2-fn-+2力,)112z、叫i+l=Pn+P)cn+i=19y-y-2+3A(-/；,_i+2fn+f (xn+i,mn+)89,、yn+l=Q+l-P+l)2.算法程序实现在 MATLAB中编程实现的Hamming算法函数为：MyHamming。功 能：以 Hamming算法求解常微分方程。调用格式：x,y=MyHammingG x0,xt,jO,N,KC)O其 中，/为 函 数 名；xO为自变量区间初值；W为自变量区间终值；298 第12章常微分方程求解yO为函数在xO处的值；N 为自变量在 xO，M 上所
34、取的点数；KC为确定是否使用改正公式；x 为输出自变量区间上所取点的x 值；y 为对应点上的函数值。Hamming算法的MATLAB程序代码如下：function x,y=MyHamming(f,xO,xt,yO,Nz KC)%用hamming方法解向量微分方程yz(t)=f(t,y(t)软函数f(x,y)：f告 自变量的初值和终值：xO,xt%函数在xO处的 值：yO出自变量在 xO,xt上取的点数：N为选择是否使用改正公式：KC生x:所取的点的x值%y:对应点上的函数值if nargin 6,KC=1;endif nargin 5 I N=0N=100;endif nargin 4yO=0
35、;endyO=yO(:)1;h=(xt-xO)/(N-l);xtO=x0+2*h;x,y=MyRunge_Kutta(fz xO,xtO,yOz 3)x=x(1:3)*x(4):h:xt 1;for k=2:4F(k-1,:)=feval(f,x(k)z y(k,:);endP=y(4,:);c=y(4,:);h34=h/3*4;KC11=KC*112/121;KC91=KC*9/121;h312=3*h*-1 2 1;带默认使用更正公式与最大迭代数默认为100考初值默认为0%使y O为行向量超步长务用Runge-Kutta得初始三点值务重新调整x向 量，使其前三点与Runge-Kutta法计
36、算的三点相同告计算f2 f3 f4%预估量初值超校正量初值称预估公式中系数务更正公式中系数与最后计算y值的公式中更正项的系数当 校正公式中各f项系数词 =0 的 y值专前推方程迭代初始化电用前推方程迭代计算n exl215y=Columns 1 through 17000000000035444 4 4Columns 18 through 3144444444444444图 1212所示作出了坐标-10,20上的卜值。此程序仅计算 10到 20 之间的p 值，依此0,2图 1 2-1 2 递推法所得y 值12.6.3用z反变换求解利用Z反变换求解差分方程多用于手工求解系统脉冲响应的过程中。MA
37、TLAB提供了符号运算工具箱，可方便地进行Z 变换和Z 反变换，进行Z 反变换的函数是iztrans,iztrans有多种调用格式，其常用的调用格式为：/=iztrans(7?)310 第12章常微分方程求解表示求尸的z 变换原函数人n)。此处仅通过一例简单介绍一下该方法。【例 12-16 利用z 变换求解差分方程实例。用单边z 反变换求解差分方程的脉冲响应(求”.0 的 值)：其中 y(-i)=。解：可知单位脉冲响应()为X()=3(N)时的输出y()。对差分方程两边进行单边变换 可 得:Y(z)-z-lY(z)-y(-1)=X(z),当 x()=b()时，X(z)=l,Y(z)=H(z)o
38、 再由y(-l)=0,故可得：用 MATLAB中的iztrans求其反变换，可得力()的关于变量的表达式。具体的程序如下所示，在命令窗口输入：syms zhz=l/(l-zA(-1);hn=iztrans(hz)输出出结果为：hn=112.7小结本章讲述了 MATLAB的求解函数和求解器及其常见用法，还讲述了采用MATLAB解决常微分方程初值问题的几种方法，包括欧拉法、预估-校正法和龙格-库塔法，最后讲述了常系数差分方程的三种基本解法。本章所举的例子很多并不是单纯用一种方法求解的，很多例子的程序中都涉及了几种算法的比较，包括计算精度及计算时间的比较。读者在选择算法解微分方程时也要考虑到这两个方面。差分方程求解法的实际应用中,filter函数多被用来求系统脉冲响应，或者作为求滤波器输出。由于现实生活中的系统为因果系统，其满足式次)=0、故用递归法不需附加条件便可求脉冲响应，其应用也很普遍。z 反变换求解多用在手工计算中。311