實驗二 MATLAB數(shù)值計算：二階電路的時域分析.doc

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