连续系统时域数值解法毕业论文

1、 题 目 基于Mathematica函数NDSolve连续系统 时域数值解法程序设计 学生姓名 刘江 学号 1110064086 所在学院 物理与电信工程学院 专业班级 电子信息科学与技术1103班 指导教师 龙姝明 完成地点 博远楼A1109 2015年5月17日 陕西理工学院毕业论文 基于Mathematica函数NDSolve连续系统时域数值解法程序设计 刘江(陕西理工学院物电学院电子信息科学与技术专业2011级 陕西 汉中 723000) 指导教师:龙姝明 摘要 采用Mathematica软件的NDSlove求连续系统的数值解，并编写求连续系统数值解的实用程序，从实例展示程序的用法和高

2、效率。 关键词 Mathematica程序; 数值解; 时域 The program design of the continuous systems time domain numerical solution Liu Jiang(Grade11,Class3,Major Electronic Information Science and TechnologyDepartment of Physics,Shannxi University of Technology,Hanzhong 723000，Shaanxi) Tutor：Long Shuming Abstract: Use the

3、NDSlove of Mathematica software to solve continuous systems numerical solution.And compile the utility program of solving continuous systems numerical solution.From the living example to diaplay programs usage and high efficency Keywords: Mathematica program;numerical solution; time domain 陕西理工学院毕业论

4、文 目录引言11 描述LTI连续系统的方法11.1 微分方程与初值的方法描述21.2 LTI连续系统解析解求解思路32 连续系统数值解的Mathematica程序设计思路52.1 系统的描述62.2 调用Mathematica的函数NDSolve求系统的数值解(响应)72.3 系统数值响应的描述73 NDslove方法求解连续系统的优缺点84 连续系统数值解的可实用Mathematica程序94.1 程序设计思路94.2 程序清单95 动态电路数值解程序应用示例105.1 时域动态电路105.1 三阶动态电路105.2 程序展示10结语11 陕西理工学院毕业论文 引言 连续系统的各种解析解法虽

5、然便于理论分析系统响应的变化趋势的系统特性，但实际系统总是多输入多输出的高阶系统，它们的解析形式的响应求解极为困难、即便较低阶系统的解析响应能够得到，其函数表示也比较复杂。而连续系统的区间数值解法本质上用的是迭代解法 ，总是能够方便、快速地的得到，之后如果企图观察其响应随时间演化的趋势，可用数值解画出其波形来观察，甚至必要时做数据拟合寻找区间解的拟合函数也是有可能的，而且数值解法还可以求一定区间上的非线性问题。1 LTI连续系统解法基本思路1.1 LTI连续系统的描述连续LTI系统的分析求解是各类电子设备、控制系统设计的基础，电子设备、控制系统设计过程中要不断地求解各类连续LTI系统在给定激励

6、(输入信号)下的响应。且众所周知，连续系统的解析解在低阶系统时比较易求，但高阶系统解析形式的响应求解极为困难，即便较低阶系统的解析响应能够得到，其函数表示也比较复杂。所以在实际中我们总是尽可能找寻系统的数值解，而连续系统的区间数值解法本质上用的是迭代解法，总是能够方便，快速的得到，之后如果企图观察其响应随时间演化的趋势，可画出其波形来观察，甚至必要时做数据拟合寻找区间解的拟合函数得到系统的演化规律也是可能的，而且数值解法还可以求解非线性问题。连续LTI可以用微分方程描述，也可以用框图或流图描述，还可以用s域系统函数（传递函数）、状态空间矩阵描述。我们的设计是编写求系统数值解的实用程序，因而重点

7、关注系统的微分方程描述方法。简单系统建立微分方程比较容易，复杂系统建立微分非常困难。为此需要将系统映射到s域，列写系统响应像函数满足的代数方程并解之，再解出s域系统函数(即传递函数)，据此可以快速写出系统的时域微分方程。1.2 LTI连续系统解析解求解思路 在实际工程系统运用中，我们总是经常需要求解n阶LTI系统的初值问题 (1-1)求解这一问题的步骤是：（1） 求解系统的零输入响应，它可由 （1-2）得到，要求解这一问题，需要先解决一元n次代数方程 （1-3）的n个根，假设它有n个单根，则有 （1-4）为了确定其中的n个迭加常数，需要将（1-2）中的常数代入（1-4）中，得到下面的n元一次代

8、数方程组 （1-5）从上面可以看出，如果n=3,这两个方程手工解都将很复杂，因此这时我们就应该寻求用编程软件Mathematica、Maple 或Matlab等软件编程或调用内部函数来解决。（2） 求问题（1-1）中的单位冲激响应 (1-6)问题（1-6）求解过程（1-2）求解过程基本一致，但可能条件有所差别，分析有 （1-7）上式中的是单位阶跃信号函数。（1-7）中的迭加系数满足n元一次代数方程组 （1-8）求h（t）的第二步求出后，由（1-1）中右端的组合系数分别乘以h(t)的各阶导数构造h（t）,即 （1-9）（3） 由h（t）和输入信号f（t）求系统的零状态响应问题（1-1）的零状态响

9、应 （1-10）（4） 求系统的完全响应 （1-11）例1 求2阶LTI连续系统的初值问题解：（1）由方程y(t)+2y(t)+5y(t)=0, y(0_)=5, y (0_)=1给出（2） 单位冲激响应h1(t)由方程上面的时域卷积计算过于复杂，可采样Laplace变换，在像域，各自的像乘积后取逆变换来计算（4）求系统的全响应 例2 求解三阶LTI连续系统初值问题解： 由于特征方程手工无法快速求解所以我们试图调用Mathematica软件的内部函数DSolve来求解这一问题，用到的程序代码为Cleary,t,sol;sol=DSolveyt+2.7yt+5.123yt+12.78yt=E(-

10、1.78t) UnitStept,y0=0,y0=0,y0=0,yt,t;yt_=yt/.sol1Plotyt,t,0,50,PlotRange-All结果计算机系统无法求解。为了能够正确求解，我们又一次试图使用NDSolve函数进行求解，其Mathematica程序代码如下:Cleary,t,sol;sol=NDSolveyt+2.7yt+5.123yt+12.78yt=E(-1.78t) UnitStept,y0=0,y0=0,y0=0,y, t,0,50;yt_=yt/.sol1Plotyt,t,0,50,PlotRange-All,AxesLabel-t,y(t)运行后输出解函数波形图

11、为图1 输出函数波形由以上求解步骤和实例可见，求解高阶LTI系统解析解必须解决两个初等代数问题:(1) 求一元n次代数方程（1-3）的根，n2时，这个问题手工很难快速求解。(2) 求两个n元一次代数方程组（1-5）和（1-8），n较大时，手工求解很困难。由此可见，对于n2的高阶LTI系统初值问题, 手工求解析解很困难。编程上机求解实践发现，调用Mathematica的内部函数求解n2的高阶LTI系统初值问题解析解也很困难，或者可以求得解析解，但求解结果是非常复杂的含复数的一大堆函数组合表示，有时候表达式长到写满一页或多页A4纸的空间。也就是说，求解高阶LTI系统初值问题解析解无论是手工求解或者

12、是编程上机求解都是令人头疼的事情。离散LTI系统初值问题的完全响应可以很方便的由迭代法求得（当然数据量很大时，肯定需要编程上机迭代求解，但算法很简单，程序代码也很简短）。如果能将高阶LTI连续系统初值问题离散化为差分方程初值问题，则由于差分方程很容易求解，如果采样间隔尽可能的小，我们可用差分方程初值问题的解，即原连续LTI系统的数值解来逼近系统的解析解。这就是我们为什么要讨论求高阶LTI连续系统数值解，并编写求数值解的实用程序设计的目的。问题已经很清楚，可以用数值解很好的逼近系统的解析解，并且工程应用和科学研究中更多用到的是系统的数值解。 线性时不变（LTI）连续系统的时域分析方法，即对于给定

13、的激励，根据描述系统的响应与系统之间关系的微分方程求得其响应的方法，其主要方法为经典解。但是利用经典解在求解微分方程的基础上讨论其零输入、零状态和全响应比较复杂，高阶解更是困难，所以求解线性时不变（LTI）连续系统就是将其映射到S域。用拉普拉斯变换法将连续系统映射到S域，写出系统函数H（s），再根据其初值和输入得出零输入、零状态和全响应。 设LTI系统的f(t),响应y（t），描述n阶系统的微分方程的一般形式可写为 （1.1-1） 式中，系数均为实数，设系统的初始状态为, 。令，。根据时域微分定理，y（t）及其各阶导数的拉普拉斯变换为 （） （1.1-2）如果f(t)是t=0时接入的，则在时f

14、(t)及其各阶导数均为零，即。因而f(t)及其各阶导数的拉普拉斯变换为 （1.1-3）取式（1.1-1）的拉普拉斯变换并将式（1.1-2）、（1.1-3）代入得 即 （1.1-4）由上式可解得 （1.1-5）式中，是方程（1.1-1）的特征多项式；，多项式和的系数仅与微分方程的系数、有关；，它也是s的多项式，其系数与和响应的各初始状态有关而与激励无关。 由式（1.1-5）可以看出，其第一项仅与初始状态有关而与输入无关，因而是零输入响应的象函数，记为；其第二项仅与激励有关而与初始状态无关，因而是零状态响应的象函数，记为。于是式（1.1-5）可写为 （1.1-6）式中，。取上式逆变换，得系统的全响

15、应 2 连续系统数值解的Mathematica程序设计思路2.1 系统的描述 对于简单系统，可以直接写出系统的微分方程。对于复杂的连续系统，直接写出微分方程比较困难，可将问题映射到s域，先导出s域系统函数H(s)，由系统函数H(s)写出微分方程,其方法是 (1) 提取H(s)的分母多项式的系数记入数组A=a0,a1,.an,提取H(s)的分子多项式的系数记入数组B=b0,b1,.bm，用到的程序语句组 B=CoefficientListNumeratorHs,s; m=LengthB-1; A=CoefficientListDenominatorHs,s; n=LengthA-1; (2) 由

16、H(s)分母分子多项式的系数写微分方程 eq=A.TableDyt,t,k,k,0,n=B.TableDft,t,k,k,0,m 例如，如图2-1所示的三阶电路系统,利用拉普拉斯变换映射到s域，方法是： 电阻不变， 图2-1 时域三阶电路 电感，电容，时域电路变成s域电路，如图2-2所示 图2-2 S域电路 由复频域电路图建立复频域代数方程，从R2两端看进去的开路像电压及等效像阻抗分别为 解得电阻R2上的像电压函数为： 代入R1=R2=1,L1=L2=1H,C=2F的值，得系统函数： 比较得到B=0.5,m=0,A=1,2,2,1,n=3。设输入为，响应为，由此列出系统的时域微分方程为 2.2

17、 调用Mathematica的函数NDSolve求系统的数值解(响应) 设微分方程已经存入变量eq中，再将初值条件存入ic变量中，例如 ic=y0=0,y0=0,y0=0,选择求解的时区为ta,tb,那么求系统数值响应的方法是 Sol=NDSolveeq,ic,y,t,ta,tb; 将插值函数形式的数值解存为一个隐函数 yt_=yt/.Sol1;至此，我们就求得了系统的数值响应(解)。例如要求图2-2电路中电阻R2两端的电压函数的数值解，所需的mathematica程序如下Cleary*,t;tf=50; eq=yt+2yt+2yt+yt=10e(-0.5t)Cos10tUnitStept;i

18、c=y0=0,y0=0,y0=0;sol=NDSolve eq,ic,y,t,ta,tb;yt_=yt/.sol12.3 系统数值响应的描述（1） 输出数值解数据表例如，如图2-2所示的电路系统，电阻R2上的电压数值解前20个数据为0.,0.,0.2,0.00959416,0.4,0.0361153,0.6,0.0388432,0.8,0.0343152,1.,0.0453352,1.2,0.0429989,1.4,0.0320751,1.6,0.034151,1.8,0.0304077,2.,0.0185072,2.2,0.0166017,2.4,0.0137429,2.6,0.003953

19、88,2.8,0.00118354,3.,0.00016305,3.2,-0.00634258,3.4,-0.00846136,3.6,-0.00768478,3.8,-0.0110255将上述数据制成如下表格：表1 数据表000.20.009594160.40.03611530.60.03884320.80.034315210.04533521.20.04299891.40.03207511.60.0341511.80.030407720.01850722.20.01660172.40.01374292.60.003953882.80.0011835430.000163053.2-0.006

20、342583.4-0.008461363.6-0.007684783.8-0.0110255我们从数值表中无法看到电阻R2上的电压随时间变化的整体趋势，为此明智的做法是将系统的数值响应表示成如图2-3所示的波形曲线，所需的程序语句为Plotyt,t,ta,tb,AxesLabel-t,y(t),PlotRange-All 图2-3 在图2-2中的电阻R2上的电压波形曲线3 NDslove方法求解连续系统的优缺点 NDslove能够极为方便的绘制函数曲线和函数曲面，求解连续系统数值解很方便，有利于做出解的图形。只要连续系统能够利用数值解求解，就一定可以用NDSlove来编写程序得出数值解波形曲线

21、。Mathematica能够求解连续系统微分方程的准确解；能求解的类型大致覆盖了人工求解的范围，功能很强，但不如人灵活（例如在隐函数和隐方程的处理方面）。输出的结果与教材上的答案可能形式不同。 另外，Mathematica只能够解数值解，不能求解连续系统的解析解。4 连续系统数值解的实用Mathematica程序4.1 程序设计思路 由电路系统的托谱结构和元件参数出发，依据基尔霍夫节点电流定律和节点电压定律，列写待求解电路元件电压或电流象函数（拉普拉斯变换）象函数满足的代数方程组，解方程组给出s域系统函数，由s域系统函数写出系统的微分方程。调用mathematica软件的NDSolve函数求系

22、统的数值解。输出数值解的数据表或者是波形曲线。4.2 程序清单对于图2-2所示电路，要求电阻R2两端电压在时区0，10上的数值解所需的mathematica程序为R1=10; R2=100; L1=0.5; L2=0.2; c=0.001;z1=R2+L2 s;z2=z1 /(z1 cs+1);z3=R1+L1 s+z2;(*U2s_=z2/z3R2/z1SubscriptU, ss;*)Hs_=z2/z3 R2/z1/Simplifynum=NumeratorHs;den=DenominatorHs;B=CoefficientListnum,s; m=LengthB-1;A=Coeffici

23、entListden,s;n=LengthA-1;ft=500Cos30t;eqL=A.TableDyt,t,j,j,0,n;eqR=B.TableDft,t,j,j,0,m;y0=-1,2,0,3;IfLengthy0n, Print初值条件数目与方程阶数数不等;ic=Table(Dyt,t,j/.t-0)=y0j+1,j,0,n-1;eq=eqL=eqR , ic/Flatten上述是由电路系统分析得出系统函数进而得出微分方程的程序，以下为求解微分方程得出数值解波形曲线的程序。 Cleary*,t; ta=0;tb=10; sol=NDSolveeq ,y,t,ta,tb; yt_=yt/

24、.sol1; Plotyt,t,ta,tb,AxesLabel-t,y(t),PlotRange-All到这里时就可以得出电路响应数值解及波形曲线了。5 动态电路数值解程序应用示例5.1 时域动态电路 图5-1为时域动态三阶电路，其中R1=10, R2=100, L1=0.5H, L2=0.2H, C=0.001F。电压为输入，电压为输出， 图5-1 三阶动态电路5.2 程序展示R1=10;R2=100;L1=0.5;L2=0.2;c=0.001;z1=R2+L2 s;z2=z1 /(z1 cs+1);z3=R1+L1 s+z2;(*U2s_=z2/z3R2/z1SubscriptU, ss;

25、*)Hs_=z2/z3 R2/z1/Simplifynum=NumeratorHsden=DenominatorHs;B=CoefficientListnum,s; m=LengthB-1;A=CoefficientListden,s;n=LengthA-1;n=Lengthden-1;ft=500Cos30t;eqL=A.TableDyt,t,j,j,0,n;eqR=B.TableDft,t,j,j,0,m;eq0=eqL=eqR;y0=-1,2,0,;ic=Table(Dyt,t,j/.t-0)=y0j+1,j,0,n-1;eq=eq0, ic/Flatten运行上述程序给出系统函数和微分

26、方程分别为：H（s）=（1106）/（1.1106+17000s+520s2+1s3） 1.1106yt+17000yt+520yt+yt=5108Cos30t,y0=-1,y0=2,y0=0调用NDSlove求解上述微分方程所需的程序 Cleary*,t; ta=0;tb=10; sol=NDSolveeq ,y,t,ta,tb; yt_=yt/.sol1; Plotyt,t,ta,tb,AxesLabel-t,y(t),PlotRange-All；画出数值解波形曲线如图5-2所示图5-2 在图5-1中电阻R2两端电压的数值解波形曲线 结语 此次毕业设计在龙姝明老师的悉心指导和严格要求下完成

27、，从最初课题的选择、方案论证和程序的编程与调试，无不凝聚着龙老师的心血和汗水，经过短暂的四年大学学习与生活我自始至终感受着老师的精心指导和无私的关怀，受益匪浅。在此向龙姝明老师表示深深的敬意和祝福。通过这次毕业设计我发现，只有理论水平提高了，才能够将课本知识与实践相整合，理论知识服务于教学实践，以增强自己的动手能力。这个设计十分有意义， 我获得很深刻的经验。通过这次毕业设计，我们知道了理论和实际的距离，也知道了理论和实际想结合的重要性，也从中得知了很多书本上无法得知的知识。我们的学习不但要立足于书本，以解决理论和实际教学中的实际问题为目的，还要以实践相结合，理论问题即实践课题，解决问题即课程研究，学生自己就是一个专家，通过自己的手来解决问题比用脑子解决问题更加深刻。学习就应该采取理论与实践结合的方式，理论的问