离散LTI系统数值解法实用程序设计初稿

离散LTI系统数值解法实用程序设计摘要：研究离散LTI系统数值解法实用程序设计，先探索连续LTI系统转化为离散LTI系统的有效方法及初值条件的转化的有效方法，再用迭代法求解指定时区上的LTI离散系统0输入响应yxi(k)和0状态响应yzs(k)，同时编写LTI连续系统转化为LTI离散系统，并用迭代方法求解系统的Mathematica程序，以LTI连续系统分析实例展示程序的高效性和可靠性，我们研究发现，数据量较大时，应用迭代编程实现更为有效更为快捷。关键词：离散系统，迭代法，mathematica程序NumericalSolutionofdiscreteLTIsystemutilitydesignedAbstract:NumericalSolutionofdiscreteLTIsystemutilitydesignedtofirstexploreeffectivemethodsofcontinuousLTIsystemintoaneffectivemethodofdiscreteLTIsystemandtheinitialconditionsoftheconversion,andthenspecifyLTIiterativemethodforsolvingdiscretesystems0inputareaontheresponseyxi(k)andzerostateresponseyzs(k),aswellaspreparefortheLTILTIcontinuoussystemintoadiscretesystemanditerativemethodsforsolvingsystemsofMathematicaprogramtoLTIcontinuoussystemsanalysisexamplesdemonstratetheefficiencyandreliabilityoftheprogram,wefoundthatlargeamountofdata,theapplicationofiterativeprogrammingmoreefficientandmoreeffectiveKeywords:Discretesystems,iterativemethod,mathematicaprogram第0页共19页引言随着数字技术以及计算机技术的飞速发展，鉴于离散系统在精度、抗干扰能力和可集成化等方面，比连续系统具有更大的优越性，因此原来对连续信号和系统的研究问题，越来越多地转化为对离散信号和系统的处理。通信和计算机设备等数字化的高科技产品渗透于人们的生活、学习、工作等诸多方面，这样，对于离散系统的分析、研究、改进成为了必不可少的课题。离散系统的响应问题是求解及分析离散系统的基础理论问题，是我们深入分析线性时不变离散系统的基础。离散LTI系统的求解方法有许多，本课题所研究的是用迭代法求解数值解。利用迭代的方法分析不借助任何变换而直接求解，直观准确。根据差分方程，用迭代的方法解出零输入响应yzi(k)，零状态响应yzs(k)。这种方法是用逐次带入来求解的，方法概念清楚，简单，对于低阶的系统手工操作就可以解出，但对于高阶系统，计算量比较大，利用计算机运算速度快、适合做重复性操作的特点，使用Mathematica软件编程实现这一过程，则更方便快捷。作为理论上的研究，此课题虽然简单，但其在基础上的意义和用途确实很不不错，为进一步深入研究奠定了很好的基础。例如在通信、计算机、自动化等工程很多领域方面都离不开对各类离散系统的分析处理，其中必定涉及高阶系统的实例。在未来的“数字化”工业发展进程中，此课题的研究方法将有更加广泛和深入的应用。第1页共19页1连续LTI系统转化为离散LTI系统的方法1.1连续系统的解析解当用数字计算机求解LTI连续系统的解析解时，或直接在系统中采用数字计算机进行求解时，对于连续低阶系统，可以通过Mathematica软件编程来实现，例如求解连续低阶系统yt+5yt+6yt=-10Cos20t,y01,y0-0.5，用Mathematica软件编程，结果如下：Cleary;eq=yt+5yt+6yt=10Cos20t,y01,y0-0.5;sol=DSolveeq,yt,t1/Expand;yt_=yt/.sol;Plotyt,t,0,5123450.20.40.60.81.0对于LTI连续低阶系统，用Mathematica软件编程可以求解它的解析解，没有多大困难。但是对于连续高阶系统，求解析解是相当困难的，例如求解yt+150000yt-12yt-5.67yt+123yt100Sin15.7t,y00,y0-4,y0-1.86,y012，用Mathematica软件编程如下：eq=yt+150000yt-12yt-5.67yt+123yt100Sin15.7t,y00,y0-4,y0-1.86,y012;sol=DSolveeq,yt,t1/Expand;yt_=yt/.sol;Plotyt,t,0,50;用这种方法求解，在Mathematica软件解的过程中出现了问题，如下面所示：RowReduce:luc：病态矩阵1.+0.,0.+0.,1.+0.,1.+0.,1.722710-7+1.6311910-21,2,1,4构成的RowReduce的结果可能包含明显的数值错误.(-3.5556410-15-3.9283710-35)-150000.t-(56.4103-1.5275210-17)-0.093707t+(56.4103-1.5276810-17)0.0468935tCos0.0809425t-(0.00161326+0.)Cos15.7t+(0.00161344_+1.6311910-21)Cos0.0809425t2Cos15.7t+(1.1991810-23-2.7377410-41)Cos0.161885tCos15.7t-(147.405_-2.6533610-17)0.0468935tSin0.0809425t+第2页共19页(4.3368110-19+2.8332410-21)Cos0.0809425tCos15.7tSin0.0809425t+(0.00161344_+0.)Cos15.7tSin0.0809425t2-(9.7419310-20+0.)Cos15.7tSin0.161885t+(9.6289210-6+3.6734210-40)Sin15.7t-(9.628910-6+9.5523410-24)Cos0.0809425t2Sin15.7t+(1.0044310-21-2.2931310-39)Cos0.161885tSin15.7t-(3.3881310-21-1.6555810-23)Cos0.0809425tSin0.0809425tSin15.7t-(9.628910-6+1.1020310-39)Sin0.0809425t2Sin15.7t+(5.815810-22+4.0783210-56)Sin0.161885tSin15.7t通过上面例子已经了解到要求解连续高阶系统的解析解是比较困难的，所以需要将它化为离散系统，来求解它在指定区间上的数值解。这样既简便又可靠，而且不会出现上面的错误。1.离散系统的数值解yt+150000yt-12yt-5.67yt+123yt100Sin15.7t,y00,y0-4,y0-1.86,y012同样地，对于上面要求解的方程,将其化为离散系统，用Mathematica对离散高阶系统进行编程，求解它在一个指定时间段上的数值解，则很容易进行求解，而且不会出现太大的错误。编程如下:Cleary;eq=yt+150000yt-12yt-5.67yt+123yt100Sin15.7t,y00,y0-4,y0-1.86,y012;sol=NDSolveeq,yt,t,0,50,MaxSteps10000001/Expand;yt_=yt/.solPlotyt,t,0,50所得到的结果如下：InterpolatingFunction0.,50.,t102030405060040020002004006008001.连续LTI系统转化为离散LTI系统由于用离散系统求解数值解更为简便、快捷和准确，故通常我们先将连续系统转化为离散系统，离散化就是导出能在采样时刻上与连续系统状态等价的离散状态方程。连续第3页共19页LTI系统转化为离散LTI系统,包含自变量的离散化、导函数的离散化和方程的离散化。例如对于方程yt+ayt+byt=cft,y0_=y0;y0_=yp0,首先是自变量的离散化，用Mathematica进行编程，如下：ta=0;tb=10;n=100;Ts=(tb-ta)/n;ts=Range0,n*Ts;然后是导函数的离散化，用Mathematica进行编程，如下：y(k)(y(k)-y(k-1)/Ts;y(k)(y(k)-2y(k-1)+y(k-2)/Ts2;y(k)(y(k)-3y(k-1)+3y(k-2)-y(k-3)/Ts3;y(k)1/Ts4(y(k)-4y(k-1)+6y(k-2)-4y(k-3)+y(k-4);.;再是方程的离散化，对于对于方程yt+ayt+byt=cft，将导函数离散化的结果代入此方程，再整理便可得到离散化的方程，用Mathematica进行编程，如下：y(k)-2y(k-1)+y(k-2)+aTs(y(k)-y(k-1)+bTs2ykcTs2fk;(1+a.Ts+b.Ts2)yk+(-2-a.Ts)yk-1+yk-2=cTs2fk;q=1/(1+a.Ts+b.Ts2);yk+(-2-aTs)qyk-1+qyk-2=cTs2qfk;p=(-2-aTs)q;=cTs2q;yk+pyk-1+qyk-2=fk;第4页共19页连续LTI系统初值条件转化为LTI离散系统初值条件的方法探索连续LTI系统初值条件转化为LTI离散系统初值条件的有效方法，首先我们是对上面导函数离散化的方程进行变形，然后赋值得到的，具体如下：对y(k)(y(k)-y(k-1)/Ts，进行变形，得到Tsy(k)y(k)-y(k-1)，然后给k赋值为-1，再移项就可得到y(-1)=y0_，对y(k)(y(k)-2y(k-1)+y(k-2)/Ts2进行变形，得到Ts2y(k)y(k)-2y(k-1)+y(k-2)，然后给k赋值为-1，再移项就可得到y(-)=y(-)y(-)+Ts2y0_，依次类推，得到如下结果：y(-1)=y0_;y(-2)=y(-1)-Tsy0_;y(-)=y(-)y(-)+Ts2y0_;.;第5页共19页3迭代法求解指定时区上的LTI离散系统零输入响应和零状态响应的编程技巧.迭代法求解指定时区上零输入响应在没有外加激励时，仅有t=0时刻的非零初始状态引起的响应。取决于初始状态和电路特性，这种相应随时间按指数规律衰减。同样对于方程yt+260300yt-12yt-6yt+123yt120Sin17t,y00,y0-4,y0-1.86,y012，用迭代法求解它在0-50s的零输入响应。即就是输入信号为的情况下求解，则上面的微分方程变为yt+260300yt-12yt-6yt+123yt，用Mathematica软件进行编程求解，如下：Cleary,u;ta=0;tb=50;n=400;Ts=(tb-ta)/n;ts=Range0,n*Ts;fn=0;q=1/(1+260300*Ts-12*Ts2-6*Ts3+123Ts4);=q*(-4-3*260300*Ts+2*12*Ts2+6*Ts3);=q*(6+3*260300*Ts-12*Ts2);=q*(-4-260300*Ts);y-1=0;y-2=4Ts;y-3=-1.86Ts2+8Ts;y-4=-3*1.86Ts2+12Ts-12Ts3;u=ConstantArray0,n+1+4;u1;4=y-4,y-3,y-2,y-1;Forj=5,jn+5,j+,uj=-*uj-1-*uj-2-*uj-3-q*uj-4;data=Transposets,u5;-1;ListPlotdata,JoinedTrue.迭代法求解指定时区上零状态响应零状态响应就是电路的储能元器件（电容、电感类元件）无初始储能，仅由外部激励作用而产生的响应。同样对于方程yt+260300yt-12yt-6yt+123yt120Sin17t,y00,y0-4,y0-1.86,y012，用迭代法求解它在0-100s的零状态响应。即就是初始状态都为的情况下求解，则y00,y0,y0,y0，用Mathematica软件进行编程求解，如下：Cleary,u;ta=0;tb=100;n=200;Ts=(tb-ta)/n;ts=Range0,n*Ts;fn=100*Sin15.7ts;第6页共19页q=1/(1+260300*Ts-12*Ts2-6*Ts3+123Ts4);=q*(-4-3*260300*Ts+2*12*Ts2+6*Ts3);=q*(6+3*260300*Ts-12*Ts2);=q*(-4-260300*Ts);y-1=0;y-2=0;y-3=0;y-4=0;u=ConstantArray0,n+1+4;u1;4=y-4,y-3,y-2,y-1;Forj=5,jn+5,j+,uj=-*uj-1-*uj-2-*uj-3-q*uj-4+q*Ts4fnj-4;data=Transposets,u5;-1;ListPlotdata,JoinedTrue20406080100.200.150.100.050.05.3迭代法求解指定时区上全响应全响应就是线性系统或电路在激励作用下产生的零状态响应与零输入响应之和。它是系统或电路在输入和初始条件共同作用下的响应。是零输入响应和零状态响应叠加的结果,也体现了线性电路的叠加性。对于方程yt+260300yt-12yt-6yt+123yt120Sin17t,y00,y0-4,y0-1.86,y012，用迭代法求解它在0-50s的全响应。即就是初始状态和输入信号都不为的情况下求解，用Mathematica软件进行编程求解，如下：Cleary,u;ta=0;tb=50;n=400;Ts=(tb-ta)/n;ts=Range0,n*Ts;fn=100*Sin15.7ts;q=1/(1+260300*Ts-12*Ts2-6*Ts3+123Ts4);=q*(-4-3*260300*Ts+2*12*Ts2+6*Ts3);=q*(6+3*260300*Ts-12*Ts2);=q*(-4-260300*Ts);y-1=0;y-2=4Ts;第7页共19页y-3=-1.86Ts2+8Ts;y-4=-3*1.86Ts2+12Ts-12Ts3;u=ConstantArray0,n+1+4;u1;4=y-4,y-3,y-2,y-1;Forj=5,jn+5,j+,uj=-*uj-1-*uj-2-*uj-3-q*uj-4+q*Ts4fnj-4;data=Transposets,u5;-1;ListPlotdata,JoinedTrue10203040506040200第8页共19页4LTI连续系统转化为LTI离散系统，迭代方法求解系统的Mathematica程序对于方程yt+150000yt-12yt-5.67yt+123yt100Sin15.7t,y00,y0-4,y0-1.86,y012，将它化为离散系统，再用迭代法求解它在0-50s的解。用Mathematica进行编程，如下：Cleary,u;ta=0;tb=50;n=200;Ts=(tb-ta)/n;ts=Range0,n*Ts;fn=100*Sin15.7ts;q=1/(1+150000*Ts-12*Ts2-5.67*Ts3+123Ts4);=q*(-4-3*150000*Ts+2*12*Ts2+5.67*Ts3);=q*(6+3*150000*Ts-12*Ts2);=q*(-4-150000*Ts);y-1=0;y-2=4Ts;y-3=-1.86Ts2+8Ts;y-4=-3*1.86Ts2+12Ts-12Ts3;u=ConstantArray0,n+1+4;u1;4=y-4,y-3,y-2,y-1;Forj=5,jn+5,j+,uj=-*uj-1-*uj-2-*uj-3-q*uj-4+q*Ts4fnj-4;在离散化的过程中，将导函数的离散化结果带入到微分方程中，并合并相同幂数的系数，将他们分别用字母表示，这样就得到离散化的差分方程为：yk+yk-1+yk-2+yk-3+qyk-4=qTs4fnk再通过Mathematica进行编程，画出它的图形：data=Transposets,u5;-1;ListPlotdata,JoinedTrue第9页共19页10203040506004002000200400600800第10页共19页5以LTI连续系统分析实例展示程序的高效性和可靠性以RLC三阶电路电路为例，具有电阻电感电容的无源二端网络如图1所示，其中：R=1，L=1H,C=2F。电压为输入，电压=150Sin18t为输出，分析该RLC2()ut三阶系统的解。图4.1RLC时域电路如图.1所示，设第一个网孔电流为，第二个网孔电流为；电容两端的电压上正1i2i下负设为，根据KVL定律，所列方程为：()cut112cccRiLuttiit21222+LRutLCiCiiiiRL整理得微分方程：22221+iiiiutLC所以输出响应为：()ut2()ti图4.2RLC频域电路利用拉普拉斯变换进行分析：电阻不变，电感，电容，时域电路Ls1Cs变成S域电路，即如图.2所示，由复频域电路图建立复频域代数方程:1OCsUUsLR其等效阻抗：0()1sCZL第11页共19页输出像函数：201()11OCRUssZULssRsCsL代入R,L,C的值，最后的输出象函数：2132UsUss由输出象函数得系统函数：321Hss根据系统函数可以列出微分方程：2yt+4yt+4yt+2yt=ft，然后通过Mathematica进行编程离散化后求解，先求零输入响应：ta=0;tb=20;n=200;Ts=(tb-ta)/n;ts=Range0,n*Ts;y(k)(y(k)-y(k-1)/Ts;y(k)(y(k)-2y(k-1)+y(k-2)/Ts2;y(k)(y(k)-3y(k-1)+3y(k-2)-y(k-3)/Ts3;q=1/(4*Ts-4*Ts2+2*Ts3+2);=q*(-6+8*Ts+4*Ts2);=q*(6+4*Ts);=q*(-2);yk+yk-1+yk-2+yk-3=qTs3fnk;y-1=0;y-2=;y-3=;u=ConstantArray0,n+1+4;u1;4=y-3,y-2,y-1;Forj=5,jn+5,j+,第12页共19页uj=-*uj-1-*uj-2-*uj-3-q*uj-4+q*Ts4fnj-4;data=Transposets,u5;-1;ListPlotdata,JoinedTrue;再用Mathematica进行编程求零输入响应：ta=0;tb=20;n=200;Ts=(tb-ta)/n;ts=Range0,n*Ts;y(k)(y(k)-y(k-1)/Ts;y(k)(y(k)-2y(k-1)+y(k-2)/Ts2;y(k)(y(k)-3y(k-1)+3y(k-2)-y(k-3)/Ts3;q=1/(4*Ts-4*Ts2-2*Ts3);=q*(-6+8*Ts+4*Ts2);=q*(6+4*Ts);=q*(-2);yk+yk-1+yk-2+yk-3=0;y-1=0;y-2=4Ts;y-3=-1.86Ts2+8Ts;y-4=-3*1.86Ts2+12Ts-12Ts3;u=ConstantArray0,n+1+4;u1;4=y-4,y-3,y-2,y-1;Forj=5,jn+5,j+,uj=-*uj-1-*uj-2-*uj-3-q*uj-4+q*Ts4fnj-4;data=Transposets,u5;-1;ListPlotdata,JoinedTrue;最后用Mathematica进行编程全响应就是零输入响应加上零状态响应：Cleary,u;ta=0;tb=50;n=400;Ts=(tb-ta)/n;第13页共19页ts=Range0,n*Ts;fn=100*Sin15.7ts;q=1/(1+260300*Ts-12*Ts2-6*Ts3+123Ts4);=q*(-4-3*260300*Ts+2*12*Ts2+6*Ts3);=q*(6+3*260300*Ts-12*Ts2);=q*(-4-260300*Ts);y-1=0;y-2=4Ts;y-3=-1.86Ts2+8Ts;y-4=-3*1.86Ts2+12Ts-12Ts3;u=ConstantArray0,n+1+4;u1;4=y-4,y-3,y-2,y-1;Forj=5,jn+5,j+,uj=-*uj-1-*