连续系统时域数值解法程序设计

连续系统时域数值解法程序设计摘要采用Mathematica软件的NDSlove求连续系统的数值解，并编写求连续系统数值解的实用程序，从实例展示程序的用法和高效率。关键词Mathematica程序;数值解;时域TheprogramdesignofthecontinuoussystemstimedomainnumericalsolutionAbstract:UsetheNDSloveofMathematicasoftwaretosolvecontinuoussystemsnumericalsolution.Andcompiletheutilityprogramofsolvingcontinuoussystemsnumericalsolution.FromthelivingexampletodiaplayprogramsusageandhighefficencyKeywords:Mathematicaprogram;numericalsolution;timedomain目录引言.11描述LTI连续系统的方法.11.1微分方程与初值的方法描述.11.2系统函数.22连续系统数值解的Mathematica程序设计思路.33复杂电路系统数值解的求解程序设计思路.43.1RLC三阶电路.43.2数学模型的建立.53.2.1列出时域连续系统电路的微分方程.53.2.2连续系统电路的时域映射到复频域（S域）.53.3NDslove方法求解连续系统的优缺点.64连续系统数值解的可实用Mathematica程序.64.1程序设计思路.64.2程序清单.65动态电路数值解程序应用示例.75.1时域动态电路.75.2程序应用展示.76结语.9参考文献.9第0页共9页引言连续系统的各种解析解法虽然便于理论分析系统响应的变化趋势的系统特性，但实际系统总是多输入多输出的高阶系统，它们的解析形式的响应求解极为困难、即便较低阶系统的解析响应能够得到，其函数表示也比较复杂1。而连续系统的区间数值解法本质上用的是迭代解法，总是能够方便、快速地的得到，之后如果企图观察其响应随时间演化的趋势，可用数值解画出其波形来观察，甚至必要时做数据拟合寻找区间解的拟合函数也是有可能的，而且数值解法还可以求一定区间上的非线性问题。第1页共9页1描述LTI连续系统的方法1.1微分方程与初值的方法描述线性时不变（LTI）连续系统的时域分析方法，即对于给定的激励，根据描述系统的响应与系统之间关系的微分方程求得其响应的方法，其主要方法为经典解。但是利用经典解在求解微分方程的基础上讨论其零输入、零状态和全响应比较复杂，高阶解更是困难，所以求解线性时不变（LTI）连续系统就是将其映射到S域2。用拉普拉斯变换法将连续系统映射到S域，写出系统函数H（s），再根据其初值和输入得出零输入、零状态和全响应3。设LTI系统的f(t),响应y（t），描述n阶系统的微分方程的一般形式可写为（1.1-1）nimjjiitfbta00)()(式中，系数均为实数，设系统的初始状态为,),.1(),.10(jbia、)0(y)(1。).1(ny令，。根据时域微分定理，y（t）及其各阶导数的拉普拉斯变换(tYsL()ftFs为（）（1.1-1()()0()iiipytsYsi=0,1n2）如果f(t)是t=0时接入的，则在时f(t)及其各阶导数均为零，即t。因而f(t)及其各阶导数的拉普拉斯变换为()0(,1)jfjm（1.1-()()jjftsFL3）取式（1.1-1）的拉普拉斯变换并将式（1.1-2）、（1.1-3）代入得1()000()()nimipjijasYsybs即（1.1-1()000()()nnimipjiijssysF4）由上式可解得（1.1-()()MsBYFsA5）式中，是方程（1.1-1）的特征多项式；，多项式和的系0()niiAsa0()mjjsb()AsB第2页共9页数仅与微分方程的系数、有关；，它也是s的多项式，其系iajb1()0()0nipiMsasy数与和响应的各初始状态有关而与激励无关。ia()0py由式（1.1-5）可以看出，其第一项仅与初始状态有关而与输入无关，因而是零输入响应的象函数，记为；其第二项仅与激励有关而与初始状态无关，因而是零状态响应()ziyt()ziYs的象函数，记为4。于是式（1.1-5）可写为zszs（1.1-()()()zizsMBsYFA6）式中，。取上式逆变换，得系统的全响应()ziMsYA()zsBFs（1.1-()()zizsyttyt7）1.2系统函数描述n阶LTI系统的微分方程一般可写为（1.2-()()00nmijijaytbft1）设是时接入的，则其零状态响应的象函数为()ft0（1.2-()zsBYFsA2）式中为激励的象函数，、分别为()Fs()ft()s（1.2-3）00()()niimjjAsasBsbs它们很容易根据微分方程写出。系统零状态响应的象函数与激励的象函数之比称为系统函数5，用表示，即()zsY()Fs()Hs第3页共9页（1.2-()zsYBHFA4）由描述系统的微分方程容易写出该系统的系统函数，反之亦然。由式（1.2-3）以及(1.2-)Hs4）可见，系统函数只与描述系统的微分方程系数、有关，即只与系统的结构、元件参()siajb数等有关，而与外界因素（激励、初始状态等）无关。引入系统函数的概念后，系统零状态响应的象函数可写为()zsyt（1.2-()zsYSF5）由所学已知，冲激响应是输入时系统的零状态响应，由于，故由ht()ft1)(tL（1.2-5）知，系统冲激响应的拉普拉斯变换()sHtL即系统的冲激响应与系统函数是拉普拉斯变换对，即()ht()（1.2-()hts6）系统的阶跃响应是输入时的零状态响应，由于，故有()gt()ftstL1)(（1.2-1()gHs7）一般情况下，若输入为，其象函数为，则零状态响应的象函数()ft()FzsYSs取上式逆变换，并由时域卷积定理，有（1.2-8)(*)(*)()()()(1111tfhsFLHLstyzzs可见时域卷积定理将连续系统的时域分析与复频域（s域）分析紧密结合起来，是系统分析方法更丰富，手段更加灵活。第4页共9页2连续系统数值解的Mathematica程序设计思路当连续系统直接写出微分方程比较困难时，即可将问题映射到S域来得到系统函数H（s），由系统函数写出微分方程比较容易，这样，对连续系统的数值解就可以通过Mathematica的NDSlove来编写程序6，从而得到数值解的波形曲线来观察。例如三阶电路系统如图2.1所示图2.1时域三阶电路通过对该三阶系统时域分析之后，得出系统微分方程为122yttytft接着便可以调用Mathematica的NDSlove如下Cleary*,t;tf=50;sol=NDSolveyt+2yt+2yt+yt=10e(-0.5t)Cos10tUnitStept,y0=0,y0=0,y0=0,yt,t,0,tf;yt_=yt/.sol1Plotyt,t,0,tf,AxesLabel-t,y(t),PlotRange-All在该程序中设好初值以后，便可以运行出数值解的波形曲线，通过观察曲线来分析得到系统的响应。第5页共9页3复杂电路系统数值解的求解程序设计思路3.1RLC三阶电路电阻、电感及电容是电路的基本元件，在交流电或电子技术中，常需要利用电阻、电感及电容元件组成不同的电路，用来改变输入正弦信号和输出正弦信号之间的相位差，可以构成各种振荡、选频电路、滤波器等。具有电阻电感电容的无源二端网络如图1所示，其中：R=1，L=1H,C=2F。现已电压为输入，电压为输出，分析RLC三阶系统。1()ut2()ut图3.1三阶时域系统图3.2S域电路3.2数学模型的建立3.2.1列出时域连续系统电路的微分方程设第一个网孔电流为，第二个网孔电流为；电容两端的电压上正下负设为1i2i()cut根据KVL定律7，所列方程为：12cccutRiLutiit计算得：22122+LRutRLCiCiii22+iR整理得微分方程：22221LiiiiutLCLC所以输出响应为：()ut2()Ruti第6页共9页3.2.2连续系统电路的时域映射到复频域（S域）利用拉普拉斯变换进行分析：电阻不变，电感，电容，时域电路变成S域电路，Ls1Cs即如图（2）所示，由复频域电路图建立复频域代数方程:1OCUssLR其等效阻抗：0()1sZLRC输出象函数：201()11OCUssULRsssLsCsR代入R,L,C的值，最后的输出象函数：2132UsUss由输出象函数得系统函数：321Hss利用系统函数列出时域的微分方程设输入为，求输出响应：1()ftut2yti12ytft通过以上对电路系统的分析，便可以调用Mathematica程序来完成电路响应数值解及波形图。3.3NDslove方法求解连续系统的优缺点NDslove能够极为方便的绘制函数曲线和函数曲面，求解连续系统数值解很方便，有利于做出解的图形。只要连续系统能够利用数值解求解，就一定可以用NDSlove来编写程序得出数值解波形曲第7页共9页线。Mathematica能够求解连续系统微分方程的准确解8；能求解的类型大致覆盖了人工求解的范围，功能很强，但不如人灵活（例如在隐函数和隐方程的处理方面）。输出的结果与教材上的答案可能形式不同。另外，Mathematica只能够解数值解，不能求解连续系统的解析解。第8页共9页4连续系统数值解的可实用Mathematica程序4.1程序设计思路之前对三阶电路系统的分析研究过程要利用程序来完成，必须调用实用的Mathematica程序。首先要明确这次解决的问题是，如何利用Mathematica程序来完成对连续系统数值解的设计，那么就要从每一步着手开始，首先分析电路系统，有哪些电子元件，明确输入与输出，然后分析电路列出所需要的参数，每一步的编程都得清楚了然9。接下来程序设计中明确怎么得到系统函数，如何来得到微分方程，以及最后要分析的数值解波形曲线。具体程序将在下面展示。4.2程序清单R1=10;R2=100;L1=0.5;L2=0.2;c=0.001;z1=R2+L2s;z2=z1/(z1cs+1);z3=R1+L1s+z2;(*U2s_=z2/z3R2/z1SubscriptU,ss;*)Hs_=z2/z3R2/z1/Simplifynum=NumeratorHsden1=DenominatorHs;den=Coefficientden1,s3,s2,s,den1/.s0/Flatten;n=Lengthden-1;ft_=500Cos30t;eq1=den.TableDyt,t,n-j,j,0,nnum*fty0=-1,2,0,3;IfLengthy0n,Print初值条件数目与方程阶数数不等;ic=Table(Dyt,t,j/.t-0)=y0j+1,j,0,n-1;eq=eq1,ic/Flatten上述是由电路系统分析得出系统函数进而得出微分方程的程序，以下为求解微分方程得出数值解波形曲线的程序。Cleary*,t;tf=50;sol=NDSolveyt+2yt+2yt+yt=10e(-0.5t)Cos10tUnitStept,y0=0,y0=0,y0=0,yt,t,0,tf;yt_=yt/.sol1Plotyt,t,0,tf,AxesLabel-t,y(t),PlotRange-All到这里时就可以得出电路响应数值解及波形曲线了。第9页共9页5动态电路数值解程序应用示例5.1时域动态电路图5.1为时域动态三阶电路7，其中R1=10,R2=100,L1=0.5H,L2=0.2H,C=0.001F。电压为输入，电压为输出，1()ut2()ut图5.1三阶动态电路图5.1三阶动态电路5.2程序展示R1=10;R2=100;L1=0.5;L2=0.2;c=0.001;z1=R2+L2s;z2=z1/(z1cs+1);z3=R1+L1s+z2;(*U2s_=z2/z3R2/z1SubscriptU,ss;*)Hs_=z2/z3R2/z1/Simplifynum=NumeratorHsden1=DenominatorHs;den=Coefficientden1,s3,s2,s,den1/.s0/Flatten;n=Lengthden-1;ft_=500Cos30t;eq1=den.TableDyt,t,n-j,j,0,nnum*fty0=-1,2,0,;ic=Table(Dyt,t,j/.t-0)=y0j+1,j,0,n-1;eq=eq1,ic/Flatten运行系统函数为：H（s）=（1106）/（1.1106+17000s+520s2+1s3）运行微分方程及初值条件如下：1.1106yt+17000yt+520yt+yt=5108Cos30t,y0=-1,y0=2,y0=0利用NDSlove求解上述微分方程10Cleary*,t;tf=10;sol=NDSolveyt+520yt+17000yt+106yt=第10页共9页5108Cos30tUnitStept,y0=-1,y0=2,y0=0,yt,t,0,tf;yt_=yt/.sol1Plotyt,t,0,tf,AxesLabel-t,y(t),PlotRange-All画出数值解波形曲线如图5.2所示246810t4020204060yt图5.2数值解波形曲线第11页共9页6结语此次毕业设计在龙姝明老师的悉心指导和严格要求下完成，从最初课题的选择、方案论证和程序的编程与调试，无不凝聚着龙老师的心血和汗水，经过短暂的四年大学学习与生活我自始至终感受着老师的精心指导和无私的关怀，受益匪浅。在此向龙姝明老师表示深深的敬意和祝福。通过这次毕业设计我发现，只有理论水平提高了，才能够将课本知识与实践相整合，理论知识服务于教学实践，以增强自己的动手能力。这个设计十分有意义，我获得很深刻的经验。通过这次毕业设计，我们知道了理论和实际的距离，也知道了理论和实际想结合的重要性，也从中得知了很多书本上无法得知的知识。我们的学习不但要立足于书本，以解决理论和实际教学中的实际问题为目的，还要以实践相结合，理论问题即实践课题，解决问题即课程研究，学生自己就是一个专家，通过自己的手来解决问题比用脑子解决问题更加深刻。学习就应该采取理论与实践结合的方式，理论的问题，也就是实践性的课题。这种做法既有助于完成理论知识的巩固，又有助于带动实践，解决实际问题，加强我们的动手能力和解决问题的能力。学如逆水行舟，不进则退。此次设计能够顺利的完成归功于老师的认真负责，当然我自己也相当努力，我更加完善和掌握了所学习的专业知识，并在设计中得以体现。正是有了他们的悉心帮助和支持，才使我的毕业论文工作顺利完成，在此向陕西理工学院，物理与电信工程学院的全体老师表示衷心的谢意。感谢他们四年来的辛勤栽培。第12页共9页参考文献1顾光旭.连续系统几种分析方法的简单讲解J.考试周刊,2010,47:206-207.2杨忠根,任蕾,陈红亮.基于因果微积分定理讨论时域分析和s域分析J.电气电子教学学报,2010,03:30-32.3赵树杰.连续线性系统时域分析中的函数匹配法J.西安电子科技大学学报,1990,01:94-101.4吴大正,杨林耀,张永瑞,等.信号与线性系统分析(第三版)M.北京：高等教育出版社,1998.191-239.5郝广波,谢里阳,李景波.连续系统可靠性分析方法研究J.机械制造,2008,09:68-71.6朝鲁,张鸿庆,唐立民.一个计算微分方程(组)对称群的Mathematica程序包及其应用J.计算物理,1997,03:12-45.7AlanV.Oppenheim.SignalsandSystems（SecondEdition）M.北京:电子工业出版社,2009.284-510.8ResearchprogressofthefractionalFouriertransforminsignal