實(shí)驗(yàn)二MATLAB數(shù)值計(jì)算：二階電路的時(shí)域分析資料

1、實(shí)驗(yàn)二 MATLAB數(shù)值計(jì)算：二階電路的時(shí)域分析、實(shí)驗(yàn)?zāi)康脑谖锢韺W(xué)和工程技術(shù)上，很多問題都可以用一個(gè)或一組常微分方程來描述，因此要解決相應(yīng)的實(shí)際問題往往需要首先求解對(duì)應(yīng)的微分方程（組）。在大多數(shù)情況下這些微分方程（組）通常是非線性的或者是超越方程（比如范德堡方程，波導(dǎo)本征值方程等），很難解析地求解（精確解），因此往往需要使用計(jì)算機(jī)數(shù)值求解（近似解）。MATLAB作為一種強(qiáng)大的科學(xué)計(jì)算語言，其在數(shù)值計(jì)算和數(shù)據(jù)的可視化方面具有無以倫比的優(yōu)勢(shì)。在解決常微分方程（組）問題上，MATLAB就提供了多種可適用于不同場(chǎng)合（如剛性和非剛性問題）下的求解器（Solver），例如ode45, ode15s，ode

2、23, ode23s等等。本次實(shí)驗(yàn)將以二階線性電路-RLC電路和二階非線性電路-范德堡電路的時(shí)域計(jì)算為例，了解和學(xué)習(xí)使用MATLAB乍為計(jì)算工具來解算復(fù)雜的微分 方程，以期達(dá)到如下幾個(gè)目的：1. 熟練使用dsolve函數(shù)解析求解常微分方程；2. 熟練運(yùn)用ode45求解器數(shù)值求解常微分方程；3. 了解狀態(tài)方程的概念，能使用MATLAB寸二階電路進(jìn)行計(jì)算和分析；二、實(shí)驗(yàn)預(yù)備知識(shí)1. 微分方程的概念未知的函數(shù)以及它的某些階的導(dǎo)數(shù)連同自變量都由一已知方程聯(lián)系在一起的方程稱為微分方程。如果未知函數(shù)是一元函數(shù)，稱為常微分方程（Ordinary differential equations ,簡(jiǎn)稱odes）

3、。n階常微分方程的一般形式（隱式）為：F（t,y,y,y,y（n） =0（ 1其中t為自變量。若方程中未知函數(shù)及其各階導(dǎo)數(shù)都是一次的，稱為線性常微分方程，否則 就是非線性微分方程，例如方程y”-二（1 - y2）y y = 0就是非線性的。2. 常微分方程的解及MATLAB旨令一階常微分方程與高階微分方程可以互化，已知一個(gè)n階常微分方程（顯式）y = f（t,y ,y , ,y ）（2）若令yy, yy, yyy（n J），可將上式化為n個(gè)一階常微分方程組:5 = fi(t, yi, y2,yn)yi = f2(t, yi, y2,yn)(3)HI HIyn = fn(t, yi, y2,yn

4、)(3)式稱為狀態(tài)方程，y!, y2,yn (即y, y； y,y(n-1)稱為狀態(tài)變量，其中yi (即y)就是常微分方程（2）式的解。（3）式中右邊的函數(shù)fi、f2、fn代表各個(gè)狀態(tài)變量的一階導(dǎo) 數(shù)的函數(shù)表達(dá)式，對(duì)于具體的方程它們有具體的形式，例如下列二階非線性微分方程：y“_J(1_y2)y y =0若令yi二y, y2 = y，可將其改寫成2個(gè)一階微分方程組(狀態(tài)方程)的形式：yi十W =(1 - y1)y2 - yi因此 fi *2； f2 -(i - y：)y2 - %。解析解只有少部分的線性常微分方程可以解析地求解(即可以算出精確的解表達(dá)式)，例如一階常系數(shù)常微分方程 dy;dt

5、= y i可以通過直接積分解出，而多數(shù)微分方程尤其是非線性 方程則很難得到解析解。有解析解的方程雖然可以手算解出，但是MATLAE也提供了 dslove指令來求方程的解析 解，其使用格式：S = dsolve (方程1方程2；,初始條件1初始條件2，自變量) 方程用字符串表示，自變量缺省值為to i階導(dǎo)數(shù)用D表示，2階導(dǎo)數(shù)用D2表示，以此類推。 S用于返回方程解析解的表達(dá)式。如果是求解方程組，則S為一個(gè)結(jié)構(gòu)體數(shù)組，它的每個(gè)域存放方程組每一個(gè)解的表達(dá)式。例i :求下列微分方程的解析解y，sin(2x) -y,y(0) =0,y(0) =i s=dsolve( D2y=si n(2*x)-y ,

6、y(0)=0 , Dy(0)=1 , x;) simplify(s)%以最簡(jiǎn)形式顯示 sans =-1/3*si n(x)*(-5+2*cos(x)% 方程的解(符號(hào)表達(dá)式)數(shù)值解對(duì)于沒有解析解的方程主要依靠計(jì)算機(jī)進(jìn)行數(shù)值求解(得到的是近似解)，例如方程y”.二(1 一 y2)y y =0就須通過計(jì)算機(jī)數(shù)值求解(結(jié)果是一系列解的數(shù)值而非表達(dá)式)。考慮n階微分方程(2)式的數(shù)值求解，它等價(jià)于一階常微分方程組(3)式。現(xiàn)將(3)式寫成矩陣形式：Y(t)二 F(t,Y(t)其中yi(t)y、Y(t)=y2(t)*=y*(5) 為n個(gè)分量的列向量(Column vector)，也稱狀態(tài)向量y(t)創(chuàng)丿

7、f2(6)n個(gè)分量的列向量，其每個(gè)兀素分別為(3)式右邊的函數(shù)表達(dá)式我們知道，微分方程要有唯一的確定解, 必須給定初值條件(to為初始時(shí)刻)：必須給定初值條件。因此方程(4)式要有確定的解r y(to)、 y(to)yi(to)Ay2(t。)丫0 =丫 (t0)=*=*yn (to)(7)Matlab提供了 ode45指令(ode是常微分方程的英文縮寫)來求解方程(4)的數(shù)值解(近似解)?；臼褂酶袷剑簍out, Yout =ode45 (odefun, tspan, Y0,options)其參數(shù)說明如下：odefun 般是用 M文件編寫的函數(shù)，odefun代表函數(shù)名，由用戶自己定義。函數(shù)返

8、回值為 式右邊的F(t,Y)= (f1, f2，,fn) T。故odefun函數(shù)的返回值應(yīng)是列向量，其最簡(jiǎn)單 的編寫格式為：其中 tfunction F = odefun (t, Y) 【作用是計(jì)算并返回4式中的F(t, Y)】時(shí)間變量，為標(biāo)量，代表計(jì)算進(jìn)程中的某時(shí)刻點(diǎn)代表狀態(tài)變量的列向量(即5式)返回值F(t, Y)，為列向量(見 6 式)ode45求解指令在計(jì)算時(shí)將會(huì)不斷地在各個(gè)時(shí)間點(diǎn)調(diào)用odefun函數(shù)，并自動(dòng)給輸入?yún)?shù)t和Y賦值。tspan指定方程的求解區(qū)間t, tf，t0是初始時(shí)刻。Y0用戶給定的初值條件，為n個(gè)分量的列向量，見(7)式。options可選項(xiàng)。一般情況下可缺省即可，若

9、用戶有特殊要求則須使用odeset指令設(shè)置options選項(xiàng)，具體用法可使用help odeset命令查詢，此處不做要求。tout列向量，輸出求解過程中區(qū)間如tf上各個(gè)計(jì)算點(diǎn)的時(shí)刻，即tout = (to t t| tf ) , to, tf上計(jì)算點(diǎn)的數(shù)目是由Matlab自動(dòng)生成。Yout輸出矩陣，其排列格式如下:Yout =ydtj%(tf)y2(t。)川 yn(Jy(to)y(to) HI(ryFo)y2(tjR川 yn (tl )I!=y(ti)卡9y(ti)m* h d(y7)qy2(tj川yn(tf)丿r0時(shí)的零輸入響應(yīng)，就必須求解(9)式。A.解析解方程(9)可以直接求解，因此有解

10、析解。下面改用MATLAB中的dsolve指令來求方程(9)的解析解，程序如下：S=dsolve(D2u=-R/L*Du-1/L/C*u,u(0)=1,Du(0)=0,t);% S為字符型數(shù)組(字串)，其值為方程9的解表達(dá)式L=0.5;C=0.02;R=12.5;% 元件參數(shù)t=0:0.01:1;%定義區(qū)間0,1上的時(shí)間序列y=eval(S);% eval串演算函數(shù)，計(jì)算字串S (即9式的解)在t時(shí)刻的值plot(t,y)%繪制電壓波形，即uc(t)的零輸入響應(yīng)B.數(shù)值解方程(9)雖然可以解析求解，但也可以使用則方程9改寫成如下狀態(tài)方程：ode45指令來近似計(jì)算。令 yi= uc, y2= d

11、uc/dt.yy2(11)Y=( yi.F =( fi，f2)T = ( y2,1/L/Uyi R/L*y2)T (12)二 y2_ 1 -LC上式也可以寫成(4)式的形式，因此狀態(tài)向量第一步：通過 M文件創(chuàng)建 ode函數(shù)，函數(shù)名 circuit_2order_odefun。functionF=circuit_2order_odefu n( t,Y)global L C R %定義全局變量 L、C、R,以實(shí)現(xiàn)參數(shù)在 MATLAB的基本工作空 %間和函數(shù)的專用空間之間數(shù)據(jù)的傳遞F=Y(2);-1/L/C*Y(1)-R/L*Y(2);% 函數(shù)返回值 F,列向量，見(12 )式注意，函數(shù)創(chuàng)建完后，必

12、須存盤，存儲(chǔ)文件名和函數(shù)名應(yīng)一致！第二步：創(chuàng)建 M腳本文件circuit_2order.mL=0.5; C=0.02; R=12.5;% 元件參數(shù)Y0=1;0;%列向量，初值條件tspan=0,1;%定義求解區(qū)間0,1tout,Yout=ode45(circuit_2order_odefu n,tspa n,Y0);plot(tout,Yout(:,1);%繪電壓波形，即uc(t)的零輸入響應(yīng)2. 二階非線性電路一范德堡(Van de Pol)電路范德堡電路由一個(gè)線性電感、一個(gè)線性電容和一個(gè)非線性電阻構(gòu)成，如圖(a)所示。非線性電阻的伏安特性如圖(b)所示，可很明顯看出電阻是非線性的。范德堡電

13、路的特性可由一個(gè)二階非線性微分方程描述(參見邱關(guān)源的電路第5版P463-P465):d2 iLdt2其中；-CT。上式就是著名的范德堡方程。該方程很難解析求解， 必須借助計(jì)算機(jī)數(shù)值求解。令力=畀2二di dt,將范德堡方程改寫成一階微分方程組（狀態(tài)方程）的形式：= y2y2(1 -y：)y2 - y(14)假定參數(shù)；=0.1，求解區(qū)間上式也可以寫成 式的形式。因此狀態(tài)向量 Y和F分別為5(0)、2（0）J丿丫0(16)75匸(f 13、Y =|，F(xiàn) = | =2也丿lf2丿嚴(yán)(1 y ) y2 - y1 丿(15)0, 100，以及初值條件請(qǐng)仿照實(shí)驗(yàn)內(nèi)容1,利用ode45指令求解范德堡方程，并

14、繪制電流Il的波形（即電流的零輸入響應(yīng)圖）。四、實(shí)驗(yàn)任務(wù)1. 輸入實(shí)驗(yàn)內(nèi)容1中提供的程序上機(jī)練習(xí)，熟悉dslove指令以及ode45指令的使用方法,并嘗試解答思考題1。2. 在任務(wù)1的基礎(chǔ)上，獨(dú)立編制實(shí)驗(yàn)內(nèi)容2的計(jì)算程序，并嘗試解答思考題3。五、思考題1. 在實(shí)驗(yàn)內(nèi)容1中，若使用ode45指令解方程9式，如何得到電流的零輸入響應(yīng)波形？2. 在實(shí)驗(yàn)內(nèi)容1中，試研究R=12 1, 3 1,，10時(shí)電壓Uc（t）的響應(yīng)波形，并把它們 繪制在一張圖上。3. 從任務(wù)2的電流波形圖（如下圖）中，你能得到什么結(jié)論？如何繪制范德堡電路的uc Il函數(shù)圖像（稱為相圖）？（提示：u _丄diL . 1 -3 _!

15、 y .1 y3 y ）UcIIIIy1 _ ys dt 3名 3附：實(shí)驗(yàn)任務(wù)2的參考圖實(shí)驗(yàn)二思考題 1函數(shù) circuit_2order_odefun.mfunction F=circuit_2order_odefun(t,Y) global L C RF=Y(2);-1/L/C*Y(1)-R/L*Y(2);腳本 circuit_2order.m global L C RL=0.5;C=0.02;R=12.5;Y0=1;0;tspan=0,1;tout,Yout=ode45(circuit_2order_odefun,tspan,Y0); plot(tout,-C*Yout(:,2);Eo

16、poll =(二)FIOAn)lo_d =0Aced42un4apolapOZJ=ooJ9 寸 po4noAno二 eo丁ued42mAOL 二& “0 0H990HJCt。J 一eqo_6 unDolos “(z) AL/iH L) A*onT(z) Al 上Mol 一eqo_6(A)un4apolapozl4=ooHd uo 口ouauruaapoapoz 4D2O報(bào)國(guó) CMend3任務(wù)二函數(shù) VandePol_3order_odefun.mfunction F=VandePol_2order_odefun(t,Y)global aF=Y(2);a*(1-Y(1)*Y(1)*Y(2)-Y(1);腳本 VandePol_3orde.mglobal a;a=01;Y0=0;1;tspan=0,100;to