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

實(shí)驗(yàn)二 MATLAB數(shù)值計(jì)算：二階電路的時(shí)域分析一、實(shí)驗(yàn)?zāi)康脑谖锢韺W(xué)和工程技術(shù)上，很多問(wèn)題都可以用一個(gè)或一組常微分方程來(lái)描述，因此要解決相應(yīng)的實(shí)際問(wèn)題往往需要首先求解對(duì)應(yīng)的微分方程(組)。在大多數(shù)情況下這些微分方程(組)通常是非線性的或者是超越方程(比如范德堡方程，波導(dǎo)本征值方程等)，很難解析地求解（精確解），因此往往需要使用計(jì)算機(jī)數(shù)值求解（近似解）。MATLAB作為一種強(qiáng)大的科學(xué)計(jì)算語(yǔ)言，其在數(shù)值計(jì)算和數(shù)據(jù)的可視化方面具有無(wú)以倫比的優(yōu)勢(shì)。在解決常微分方程(組)問(wèn)題上，MATLAB就提供了多種可適用于不同場(chǎng)合(如剛性和非剛性問(wèn)題)下的求解器(Solver)，例如ode45，ode15s，ode23，ode23s等等。本次實(shí)驗(yàn)將以二階線性電路-RLC電路和二階非線性電路-范德堡電路的時(shí)域計(jì)算為例，了解和學(xué)習(xí)使用MATLAB作為計(jì)算工具來(lái)解算復(fù)雜的微分方程，以期達(dá)到如下幾個(gè)目的：1. 熟練使用dsolve函數(shù)解析求解常微分方程；2. 熟練運(yùn)用ode45求解器數(shù)值求解常微分方程； 3. 了解狀態(tài)方程的概念，能使用MATLAB對(duì)二階電路進(jìn)行計(jì)算和分析；二、實(shí)驗(yàn)預(yù)備知識(shí)1微分方程的概念未知的函數(shù)以及它的某些階的導(dǎo)數(shù)連同自變量都由一已知方程聯(lián)系在一起的方程稱(chēng)為微分方程。如果未知函數(shù)是一元函數(shù)，稱(chēng)為常微分方程（Ordinary differential equations，簡(jiǎn)稱(chēng)odes）。n階常微分方程的一般形式（隱式）為： （1）其中t為自變量。若方程中未知函數(shù)及其各階導(dǎo)數(shù)都是一次的，稱(chēng)為線性常微分方程，否則就是非線性微分方程，例如方程就是非線性的。2常微分方程的解及MATLAB指令一階常微分方程與高階微分方程可以互化，已知一個(gè)n階常微分方程（顯式）： （2）若令，可將上式化為n個(gè)一階常微分方程組：（3）（3）式稱(chēng)為狀態(tài)方程，y1, y2, ,yn（即y, y, y, , y(n-1) ）稱(chēng)為狀態(tài)變量，其中y1（即y）就是常微分方程(2)式的解。(3)式中右邊的函數(shù)f1、f2、fn代表各個(gè)狀態(tài)變量的一階導(dǎo)數(shù)的函數(shù)表達(dá)式，對(duì)于具體的方程它們有具體的形式，例如下列二階非線性微分方程：若令，可將其改寫(xiě)成2個(gè)一階微分方程組(狀態(tài)方程)的形式：因此。l 解析解只有少部分的線性常微分方程可以解析地求解（即可以算出精確的解表達(dá)式），例如一階常系數(shù)常微分方程可以通過(guò)直接積分解出，而多數(shù)微分方程尤其是非線性方程則很難得到解析解。有解析解的方程雖然可以手算解出，但是MATLAB也提供了dslove指令來(lái)求方程的解析解，其使用格式：S = dsolve (方程1, 方程2,初始條件1，初始條件2 ,自變量) 方程用字符串表示，自變量缺省值為t。1階導(dǎo)數(shù)用D表示，2階導(dǎo)數(shù)用D2表示，以此類(lèi)推。S用于返回方程解析解的表達(dá)式。如果是求解方程組，則S為一個(gè)結(jié)構(gòu)體數(shù)組，它的每個(gè)域存放方程組每一個(gè)解的表達(dá)式。例1：求下列微分方程的解析解 s=dsolve(D2y=sin(2*x)-y,y(0)=0,Dy(0)=1,x); simplify(s) %以最簡(jiǎn)形式顯示sans =-1/3*sin(x)*(-5+2*cos(x) % 方程的解（符號(hào)表達(dá)式）l 數(shù)值解對(duì)于沒(méi)有解析解的方程主要依靠計(jì)算機(jī)進(jìn)行數(shù)值求解（得到的是近似解），例如方程就須通過(guò)計(jì)算機(jī)數(shù)值求解（結(jié)果是一系列解的數(shù)值而非表達(dá)式）?？紤]n階微分方程（2）式的數(shù)值求解，它等價(jià)于一階常微分方程組(3)式?，F(xiàn)將(3)式寫(xiě)成矩陣形式：其中 為n個(gè)分量的列向量(Column vector)，也稱(chēng)狀態(tài)向量 n個(gè)分量的列向量，其每個(gè)元素分別為(3)式右邊的函數(shù)表達(dá)式 我們知道，微分方程要有唯一的確定解，必須給定初值條件。因此方程(4)式要有確定的解必須給定初值條件(t0為初始時(shí)刻)： （7）Matlab提供了ode45指令（ode是常微分方程的英文縮寫(xiě)）來(lái)求解方程(4)的數(shù)值解（近似解）。基本使用格式：tout, Yout =ode45 (odefun, tspan, Y0,options) 其參數(shù)說(shuō)明如下： odefun 一般是用M文件編寫(xiě)的函數(shù)，odefun代表函數(shù)名，由用戶(hù)自己定義。函數(shù)返回值為(4)式右邊的F(t,Y)=（f1, f2 , fn）T。故odefun函數(shù)的返回值應(yīng)是列向量，其最簡(jiǎn)單的編寫(xiě)格式為： function F = odefun (t, Y) 作用是計(jì)算并返回4式中的F(t, Y) 】其中 t 時(shí)間變量，為標(biāo)量，代表計(jì)算進(jìn)程中的某時(shí)刻點(diǎn)Y 代表狀態(tài)變量的列向量(即5式)F 返回值F(t, Y)，為列向量（見(jiàn)6式）ode45求解指令在計(jì)算時(shí)將會(huì)不斷地在各個(gè)時(shí)間點(diǎn)調(diào)用odefun函數(shù)，并自動(dòng)給輸入?yún)?shù)t和Y賦值。 tspan 指定方程的求解區(qū)間t0, tf，t0是初始時(shí)刻。 Y0 用戶(hù)給定的初值條件，為n個(gè)分量的列向量，見(jiàn)(7)式。 options 可選項(xiàng)。一般情況下可缺省即可，若用戶(hù)有特殊要求則須使用odeset指令設(shè)置options選項(xiàng)，具體用法可使用help odeset命令查詢(xún)，此處不做要求。 tout 列向量，輸出求解過(guò)程中區(qū)間t0, tf上各個(gè)計(jì)算點(diǎn)的時(shí)刻，即, t0, tf上計(jì)算點(diǎn)的數(shù)目是由Matlab自動(dòng)生成。 Yout 輸出矩陣，其排列格式如下：Yout的第1列代表的是狀態(tài)變量y1在（依次為tout的每個(gè)元素）各個(gè)時(shí)刻的值，由于y1=y，所以Yout的第一列就是待求方程的數(shù)值解，它顯然是一系列離散的y(t)值：；將Yout第一列提取出來(lái)并用plot(tout,Yout(:,1)指令即可繪制解y(t)的函數(shù)圖像。Yout的第2列是y的一階導(dǎo)數(shù)的數(shù)值解，第3列是y的二階導(dǎo)數(shù)的數(shù)值解。三、實(shí)驗(yàn)內(nèi)容1. 二階線性電路RLC回路的零輸入響應(yīng)當(dāng)電路中含有二個(gè)動(dòng)態(tài)元件(如電感、電容)時(shí)，建立的電路方程為二階微分方程，這樣的電路稱(chēng)為二階電路，如果微分方程是線性的，則為線性電路，若為非線性方程，則是非線性電路，如范德堡電路。所謂零輸入響應(yīng)指的是電路中無(wú)外加的激勵(lì)源，僅由動(dòng)態(tài)元件初始儲(chǔ)能所產(chǎn)生的響應(yīng)。 考慮如圖所示的RLC電路，假設(shè)電容原已充電，根據(jù)電路理論，此二階電路的零輸入響應(yīng)可用如下二階線性微分方程描述：其中uc代表電容電壓。給定初值條件：元件參數(shù)L=0.5H, R=12.5W, C=0.02F。要得到t0時(shí)的零輸入響應(yīng)，就必須求解(9)式。A. 解析解方程(9)可以直接求解，因此有解析解。下面改用MATLAB中的dsolve指令來(lái)求方程(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)雖然可以解析求解，但也可以使用ode45指令來(lái)近似計(jì)算。令y1= uc, y2= duc/dt,則方程9改寫(xiě)成如下?tīng)顟B(tài)方程：上式也可以寫(xiě)成(4)式的形式，因此狀態(tài)向量 Y=（y1，y2）T， F =（f1，f2）T =（y2，1/L/C*y1R/L*y2）T （12）第一步：通過(guò)M文件創(chuàng)建ode函數(shù)，函數(shù)名circuit_2order_odefun。function F=circuit_2order_odefun(t,Y)global L C R % 定義全局變量L、C、R，以實(shí)現(xiàn)參數(shù)在MATLAB的基本工作空% 間和函數(shù)的專(zhuān)用空間之間數(shù)據(jù)的傳遞F=Y(2);-1/L/C*Y(1)-R/L*Y(2); % 函數(shù)返回值F, 列向量，見(jiàn)（12）式注意，函數(shù)創(chuàng)建完后，必須存盤(pán)，存儲(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_odefun,tspan,Y0);plot(tout,Yout(:,1); % 繪電壓波形，即uc(t)的零輸入響應(yīng)2. 二階非線性電路范德堡(Van de Pol)電路范德堡電路由一個(gè)線性電感、一個(gè)線性電容和一個(gè)非線性電阻構(gòu)成，如圖(a)所示。非線性電阻的伏安特性如圖(b)所示，可很明顯看出電阻是非線性的。范德堡電路的特性可由一個(gè)二階非線性微分方程描述（參見(jiàn)邱關(guān)源的電路第5版P463-P465）：其中。上式就是著名的范德堡方程。該方程很難解析求解，必須借助計(jì)算機(jī)數(shù)值求解。令，將范德堡方程改寫(xiě)成一階微分方程組(狀態(tài)方程)的形式： 上式也可以寫(xiě)成(4)式的形式。因此狀態(tài)向量Y和F分別為， （15）假定參數(shù)e =0.1，求解區(qū)間0, 100，以及初值條件請(qǐng)仿照實(shí)驗(yàn)內(nèi)容1，利用ode45指令求解范德堡方程，并繪制電流iL的波形（即電流的零輸入響應(yīng)圖）。四、實(shí)驗(yàn)任務(wù)1. 輸入實(shí)驗(yàn)內(nèi)容1中提供的程序上機(jī)練習(xí)，熟悉dslove指令以及ode45指令的使用方法，并嘗試解答思考題1。2. 在任務(wù)1的基礎(chǔ)上