連續(xù)系統(tǒng)時域數(shù)值解法畢業(yè)論文

1、 題 目 基于Mathematica函數(shù)NDSolve連續(xù)系統(tǒng) 時域數(shù)值解法程序設(shè)計 學(xué)生姓名 劉江 學(xué)號 1110064086 所在學(xué)院 物理與電信工程學(xué)院 專業(yè)班級 電子信息科學(xué)與技術(shù)1103班 指導(dǎo)教師 龍姝明 完成地點 博遠樓A1109 2015年5月17日 陜西理工學(xué)院畢業(yè)論文 基于Mathematica函數(shù)NDSolve連續(xù)系統(tǒng)時域數(shù)值解法程序設(shè)計 劉江(陜西理工學(xué)院物電學(xué)院電子信息科學(xué)與技術(shù)專業(yè)2011級 陜西 漢中 723000) 指導(dǎo)教師:龍姝明 摘要 采用Mathematica軟件的NDSlove求連續(xù)系統(tǒng)的數(shù)值解，并編寫求連續(xù)系統(tǒng)數(shù)值解的實用程序，從實例展示程序的用法和高

2、效率。 關(guān)鍵詞 Mathematica程序; 數(shù)值解; 時域 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 陜西理工學(xué)院畢業(yè)論

4、文 目錄引言11 描述LTI連續(xù)系統(tǒng)的方法11.1 微分方程與初值的方法描述21.2 LTI連續(xù)系統(tǒng)解析解求解思路32 連續(xù)系統(tǒng)數(shù)值解的Mathematica程序設(shè)計思路52.1 系統(tǒng)的描述62.2 調(diào)用Mathematica的函數(shù)NDSolve求系統(tǒng)的數(shù)值解(響應(yīng))72.3 系統(tǒng)數(shù)值響應(yīng)的描述73 NDslove方法求解連續(xù)系統(tǒng)的優(yōu)缺點84 連續(xù)系統(tǒng)數(shù)值解的可實用Mathematica程序94.1 程序設(shè)計思路94.2 程序清單95 動態(tài)電路數(shù)值解程序應(yīng)用示例105.1 時域動態(tài)電路105.1 三階動態(tài)電路105.2 程序展示10結(jié)語11 陜西理工學(xué)院畢業(yè)論文 引言 連續(xù)系統(tǒng)的各種解析解法雖

5、然便于理論分析系統(tǒng)響應(yīng)的變化趨勢的系統(tǒng)特性，但實際系統(tǒng)總是多輸入多輸出的高階系統(tǒng)，它們的解析形式的響應(yīng)求解極為困難、即便較低階系統(tǒng)的解析響應(yīng)能夠得到，其函數(shù)表示也比較復(fù)雜。而連續(xù)系統(tǒng)的區(qū)間數(shù)值解法本質(zhì)上用的是迭代解法 ，總是能夠方便、快速地的得到，之后如果企圖觀察其響應(yīng)隨時間演化的趨勢，可用數(shù)值解畫出其波形來觀察，甚至必要時做數(shù)據(jù)擬合尋找區(qū)間解的擬合函數(shù)也是有可能的，而且數(shù)值解法還可以求一定區(qū)間上的非線性問題。1 LTI連續(xù)系統(tǒng)解法基本思路1.1 LTI連續(xù)系統(tǒng)的描述連續(xù)LTI系統(tǒng)的分析求解是各類電子設(shè)備、控制系統(tǒng)設(shè)計的基礎(chǔ)，電子設(shè)備、控制系統(tǒng)設(shè)計過程中要不斷地求解各類連續(xù)LTI系統(tǒng)在給定激勵

6、(輸入信號)下的響應(yīng)。且眾所周知，連續(xù)系統(tǒng)的解析解在低階系統(tǒng)時比較易求，但高階系統(tǒng)解析形式的響應(yīng)求解極為困難，即便較低階系統(tǒng)的解析響應(yīng)能夠得到，其函數(shù)表示也比較復(fù)雜。所以在實際中我們總是盡可能找尋系統(tǒng)的數(shù)值解，而連續(xù)系統(tǒng)的區(qū)間數(shù)值解法本質(zhì)上用的是迭代解法，總是能夠方便，快速的得到，之后如果企圖觀察其響應(yīng)隨時間演化的趨勢，可畫出其波形來觀察，甚至必要時做數(shù)據(jù)擬合尋找區(qū)間解的擬合函數(shù)得到系統(tǒng)的演化規(guī)律也是可能的，而且數(shù)值解法還可以求解非線性問題。連續(xù)LTI可以用微分方程描述，也可以用框圖或流圖描述，還可以用s域系統(tǒng)函數(shù)（傳遞函數(shù)）、狀態(tài)空間矩陣描述。我們的設(shè)計是編寫求系統(tǒng)數(shù)值解的實用程序，因而重點

7、關(guān)注系統(tǒng)的微分方程描述方法。簡單系統(tǒng)建立微分方程比較容易，復(fù)雜系統(tǒng)建立微分非常困難。為此需要將系統(tǒng)映射到s域，列寫系統(tǒng)響應(yīng)像函數(shù)滿足的代數(shù)方程并解之，再解出s域系統(tǒng)函數(shù)(即傳遞函數(shù))，據(jù)此可以快速寫出系統(tǒng)的時域微分方程。1.2 LTI連續(xù)系統(tǒng)解析解求解思路 在實際工程系統(tǒng)運用中，我們總是經(jīng)常需要求解n階LTI系統(tǒng)的初值問題 (1-1)求解這一問題的步驟是：（1） 求解系統(tǒng)的零輸入響應(yīng)，它可由 （1-2）得到，要求解這一問題，需要先解決一元n次代數(shù)方程 （1-3）的n個根，假設(shè)它有n個單根，則有 （1-4）為了確定其中的n個迭加常數(shù)，需要將（1-2）中的常數(shù)代入（1-4）中，得到下面的n元一次代

8、數(shù)方程組 （1-5）從上面可以看出，如果n=3,這兩個方程手工解都將很復(fù)雜，因此這時我們就應(yīng)該尋求用編程軟件Mathematica、Maple 或Matlab等軟件編程或調(diào)用內(nèi)部函數(shù)來解決。（2） 求問題（1-1）中的單位沖激響應(yīng) (1-6)問題（1-6）求解過程（1-2）求解過程基本一致，但可能條件有所差別，分析有 （1-7）上式中的是單位階躍信號函數(shù)。（1-7）中的迭加系數(shù)滿足n元一次代數(shù)方程組 （1-8）求h（t）的第二步求出后，由（1-1）中右端的組合系數(shù)分別乘以h(t)的各階導(dǎo)數(shù)構(gòu)造h（t）,即 （1-9）（3） 由h（t）和輸入信號f（t）求系統(tǒng)的零狀態(tài)響應(yīng)問題（1-1）的零狀態(tài)響

9、應(yīng) （1-10）（4） 求系統(tǒng)的完全響應(yīng) （1-11）例1 求2階LTI連續(xù)系統(tǒng)的初值問題解：（1）由方程y(t)+2y(t)+5y(t)=0, y(0_)=5, y (0_)=1給出（2） 單位沖激響應(yīng)h1(t)由方程上面的時域卷積計算過于復(fù)雜，可采樣Laplace變換，在像域，各自的像乘積后取逆變換來計算（4）求系統(tǒng)的全響應(yīng) 例2 求解三階LTI連續(xù)系統(tǒng)初值問題解： 由于特征方程手工無法快速求解所以我們試圖調(diào)用Mathematica軟件的內(nèi)部函數(shù)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結(jié)果計算機系統(tǒng)無法求解。為了能夠正確求解，我們又一次試圖使用NDSolve函數(shù)進行求解，其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)運行后輸出解函數(shù)波形圖

11、為圖1 輸出函數(shù)波形由以上求解步驟和實例可見，求解高階LTI系統(tǒng)解析解必須解決兩個初等代數(shù)問題:(1) 求一元n次代數(shù)方程（1-3）的根，n2時，這個問題手工很難快速求解。(2) 求兩個n元一次代數(shù)方程組（1-5）和（1-8），n較大時，手工求解很困難。由此可見，對于n2的高階LTI系統(tǒng)初值問題, 手工求解析解很困難。編程上機求解實踐發(fā)現(xiàn)，調(diào)用Mathematica的內(nèi)部函數(shù)求解n2的高階LTI系統(tǒng)初值問題解析解也很困難，或者可以求得解析解，但求解結(jié)果是非常復(fù)雜的含復(fù)數(shù)的一大堆函數(shù)組合表示，有時候表達式長到寫滿一頁或多頁A4紙的空間。也就是說，求解高階LTI系統(tǒng)初值問題解析解無論是手工求解或者

12、是編程上機求解都是令人頭疼的事情。離散LTI系統(tǒng)初值問題的完全響應(yīng)可以很方便的由迭代法求得（當(dāng)然數(shù)據(jù)量很大時，肯定需要編程上機迭代求解，但算法很簡單，程序代碼也很簡短）。如果能將高階LTI連續(xù)系統(tǒng)初值問題離散化為差分方程初值問題，則由于差分方程很容易求解，如果采樣間隔盡可能的小，我們可用差分方程初值問題的解，即原連續(xù)LTI系統(tǒng)的數(shù)值解來逼近系統(tǒng)的解析解。這就是我們?yōu)槭裁匆懻撉蟾唠ALTI連續(xù)系統(tǒng)數(shù)值解，并編寫求數(shù)值解的實用程序設(shè)計的目的。問題已經(jīng)很清楚，可以用數(shù)值解很好的逼近系統(tǒng)的解析解，并且工程應(yīng)用和科學(xué)研究中更多用到的是系統(tǒng)的數(shù)值解。 線性時不變（LTI）連續(xù)系統(tǒng)的時域分析方法，即對于給定

13、的激勵，根據(jù)描述系統(tǒng)的響應(yīng)與系統(tǒng)之間關(guān)系的微分方程求得其響應(yīng)的方法，其主要方法為經(jīng)典解。但是利用經(jīng)典解在求解微分方程的基礎(chǔ)上討論其零輸入、零狀態(tài)和全響應(yīng)比較復(fù)雜，高階解更是困難，所以求解線性時不變（LTI）連續(xù)系統(tǒng)就是將其映射到S域。用拉普拉斯變換法將連續(xù)系統(tǒng)映射到S域，寫出系統(tǒng)函數(shù)H（s），再根據(jù)其初值和輸入得出零輸入、零狀態(tài)和全響應(yīng)。 設(shè)LTI系統(tǒng)的f(t),響應(yīng)y（t），描述n階系統(tǒng)的微分方程的一般形式可寫為 （1.1-1） 式中，系數(shù)均為實數(shù)，設(shè)系統(tǒng)的初始狀態(tài)為, 。令，。根據(jù)時域微分定理，y（t）及其各階導(dǎo)數(shù)的拉普拉斯變換為 （） （1.1-2）如果f(t)是t=0時接入的，則在時f

14、(t)及其各階導(dǎo)數(shù)均為零，即。因而f(t)及其各階導(dǎo)數(shù)的拉普拉斯變換為 （1.1-3）取式（1.1-1）的拉普拉斯變換并將式（1.1-2）、（1.1-3）代入得 即 （1.1-4）由上式可解得 （1.1-5）式中，是方程（1.1-1）的特征多項式；，多項式和的系數(shù)僅與微分方程的系數(shù)、有關(guān)；，它也是s的多項式，其系數(shù)與和響應(yīng)的各初始狀態(tài)有關(guān)而與激勵無關(guān)。 由式（1.1-5）可以看出，其第一項僅與初始狀態(tài)有關(guān)而與輸入無關(guān)，因而是零輸入響應(yīng)的象函數(shù)，記為；其第二項僅與激勵有關(guān)而與初始狀態(tài)無關(guān)，因而是零狀態(tài)響應(yīng)的象函數(shù)，記為。于是式（1.1-5）可寫為 （1.1-6）式中，。取上式逆變換，得系統(tǒng)的全響

15、應(yīng) 2 連續(xù)系統(tǒng)數(shù)值解的Mathematica程序設(shè)計思路2.1 系統(tǒng)的描述 對于簡單系統(tǒng)，可以直接寫出系統(tǒng)的微分方程。對于復(fù)雜的連續(xù)系統(tǒng)，直接寫出微分方程比較困難，可將問題映射到s域，先導(dǎo)出s域系統(tǒng)函數(shù)H(s)，由系統(tǒng)函數(shù)H(s)寫出微分方程,其方法是 (1) 提取H(s)的分母多項式的系數(shù)記入數(shù)組A=a0,a1,.an,提取H(s)的分子多項式的系數(shù)記入數(shù)組B=b0,b1,.bm，用到的程序語句組 B=CoefficientListNumeratorHs,s; m=LengthB-1; A=CoefficientListDenominatorHs,s; n=LengthA-1; (2) 由

16、H(s)分母分子多項式的系數(shù)寫微分方程 eq=A.TableDyt,t,k,k,0,n=B.TableDft,t,k,k,0,m 例如，如圖2-1所示的三階電路系統(tǒng),利用拉普拉斯變換映射到s域，方法是： 電阻不變， 圖2-1 時域三階電路 電感，電容，時域電路變成s域電路，如圖2-2所示 圖2-2 S域電路 由復(fù)頻域電路圖建立復(fù)頻域代數(shù)方程，從R2兩端看進去的開路像電壓及等效像阻抗分別為 解得電阻R2上的像電壓函數(shù)為： 代入R1=R2=1,L1=L2=1H,C=2F的值，得系統(tǒng)函數(shù)： 比較得到B=0.5,m=0,A=1,2,2,1,n=3。設(shè)輸入為，響應(yīng)為，由此列出系統(tǒng)的時域微分方程為 2.2

17、 調(diào)用Mathematica的函數(shù)NDSolve求系統(tǒng)的數(shù)值解(響應(yīng)) 設(shè)微分方程已經(jīng)存入變量eq中，再將初值條件存入ic變量中，例如 ic=y0=0,y0=0,y0=0,選擇求解的時區(qū)為ta,tb,那么求系統(tǒng)數(shù)值響應(yīng)的方法是 Sol=NDSolveeq,ic,y,t,ta,tb; 將插值函數(shù)形式的數(shù)值解存為一個隱函數(shù) yt_=yt/.Sol1;至此，我們就求得了系統(tǒng)的數(shù)值響應(yīng)(解)。例如要求圖2-2電路中電阻R2兩端的電壓函數(shù)的數(shù)值解，所需的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 系統(tǒng)數(shù)值響應(yīng)的描述（1） 輸出數(shù)值解數(shù)據(jù)表例如，如圖2-2所示的電路系統(tǒng)，電阻R2上的電壓數(shù)值解前20個數(shù)據(jù)為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將上述數(shù)據(jù)制成如下表格：表1 數(shù)據(jù)表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我們從數(shù)值表中無法看到電阻R2上的電壓隨時間變化的整體趨勢，為此明智的做法是將系統(tǒng)的數(shù)值響應(yīng)表示成如圖2-3所示的波形曲線，所需的程序語句為Plotyt,t,ta,tb,AxesLabel-t,y(t),PlotRange-All 圖2-3 在圖2-2中的電阻R2上的電壓波形曲線3 NDslove方法求解連續(xù)系統(tǒng)的優(yōu)缺點 NDslove能夠極為方便的繪制函數(shù)曲線和函數(shù)曲面，求解連續(xù)系統(tǒng)數(shù)值解很方便，有利于做出解的圖形。只要連續(xù)系統(tǒng)能夠利用數(shù)值解求解，就一定可以用NDSlove來編寫程序得出數(shù)值解波形曲線

21、。Mathematica能夠求解連續(xù)系統(tǒng)微分方程的準(zhǔn)確解；能求解的類型大致覆蓋了人工求解的范圍，功能很強，但不如人靈活（例如在隱函數(shù)和隱方程的處理方面）。輸出的結(jié)果與教材上的答案可能形式不同。 另外，Mathematica只能夠解數(shù)值解，不能求解連續(xù)系統(tǒng)的解析解。4 連續(xù)系統(tǒng)數(shù)值解的實用Mathematica程序4.1 程序設(shè)計思路 由電路系統(tǒng)的托譜結(jié)構(gòu)和元件參數(shù)出發(fā)，依據(jù)基爾霍夫節(jié)點電流定律和節(jié)點電壓定律，列寫待求解電路元件電壓或電流象函數(shù)（拉普拉斯變換）象函數(shù)滿足的代數(shù)方程組，解方程組給出s域系統(tǒng)函數(shù)，由s域系統(tǒng)函數(shù)寫出系統(tǒng)的微分方程。調(diào)用mathematica軟件的NDSolve函數(shù)求系

22、統(tǒng)的數(shù)值解。輸出數(shù)值解的數(shù)據(jù)表或者是波形曲線。4.2 程序清單對于圖2-2所示電路，要求電阻R2兩端電壓在時區(qū)0，10上的數(shù)值解所需的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初值條件數(shù)目與方程階數(shù)數(shù)不等;ic=Table(Dyt,t,j/.t-0)=y0j+1,j,0,n-1;eq=eqL=eqR , ic/Flatten上述是由電路系統(tǒng)分析得出系統(tǒng)函數(shù)進而得出微分方程的程序，以下為求解微分方程得出數(shù)值解波形曲線的程序。 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到這里時就可以得出電路響應(yīng)數(shù)值解及波形曲線了。5 動態(tài)電路數(shù)值解程序應(yīng)用示例5.1 時域動態(tài)電路 圖5-1為時域動態(tài)三階電路，其中R1=10, R2=100, L1=0.5H, L2=0.2H, C=0.001F。電壓為輸入，電壓為輸出， 圖5-1 三階動態(tài)電路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運行上述程序給出系統(tǒng)函數(shù)和微分

26、方程分別為：H（s）=（1106）/（1.1106+17000s+520s2+1s3） 1.1106yt+17000yt+520yt+yt=5108Cos30t,y0=-1,y0=2,y0=0調(diào)用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；畫出數(shù)值解波形曲線如圖5-2所示圖5-2 在圖5-1中電阻R2兩端電壓的數(shù)值解波形曲線 結(jié)語 此次畢業(yè)設(shè)計在龍姝明老師的悉心指導(dǎo)和嚴(yán)格要求下完成

27、，從最初課題的選擇、方案論證和程序的編程與調(diào)試，無不凝聚著龍老師的心血和汗水，經(jīng)過短暫的四年大學(xué)學(xué)習(xí)與生活我自始至終感受著老師的精心指導(dǎo)和無私的關(guān)懷，受益匪淺。在此向龍姝明老師表示深深的敬意和祝福。通過這次畢業(yè)設(shè)計我發(fā)現(xiàn)，只有理論水平提高了，才能夠?qū)⒄n本知識與實踐相整合，理論知識服務(wù)于教學(xué)實踐，以增強自己的動手能力。這個設(shè)計十分有意義， 我獲得很深刻的經(jīng)驗。通過這次畢業(yè)設(shè)計，我們知道了理論和實際的距離，也知道了理論和實際想結(jié)合的重要性，也從中得知了很多書本上無法得知的知識。我們的學(xué)習(xí)不但要立足于書本，以解決理論和實際教學(xué)中的實際問題為目的，還要以實踐相結(jié)合，理論問題即實踐課題，解決問題即課程研究，學(xué)生自己就是一個專家，通過自己的手來解決問題比用腦子解決問題更加深刻。學(xué)習(xí)就應(yīng)該采取理論與實踐結(jié)合的方式，理論的問