实验二MATLAB数值计算：二阶电路的时域分析资料

1、实验二 MATLAB数值计算：二阶电路的时域分析、实验目的在物理学和工程技术上，很多问题都可以用一个或一组常微分方程来描述，因此要解决相应的实际问题往往需要首先求解对应的微分方程（组）。在大多数情况下这些微分方程（组）通常是非线性的或者是超越方程（比如范德堡方程，波导本征值方程等），很难解析地求解（精确解），因此往往需要使用计算机数值求解（近似解）。MATLAB作为一种强大的科学计算语言，其在数值计算和数据的可视化方面具有无以伦比的优势。在解决常微分方程（组）问题上，MATLAB就提供了多种可适用于不同场合（如刚性和非刚性问题）下的求解器（Solver），例如ode45, ode15s，ode

2、23, ode23s等等。本次实验将以二阶线性电路-RLC电路和二阶非线性电路-范德堡电路的时域计算为例，了解和学习使用MATLAB乍为计算工具来解算复杂的微分 方程，以期达到如下几个目的：1. 熟练使用dsolve函数解析求解常微分方程；2. 熟练运用ode45求解器数值求解常微分方程；3. 了解状态方程的概念，能使用MATLAB寸二阶电路进行计算和分析；二、实验预备知识1. 微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程（Ordinary differential equations ,简称odes）

3、。n阶常微分方程的一般形式（隐式）为：F（t,y,y,y,y（n） =0（ 1其中t为自变量。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，否则 就是非线性微分方程，例如方程y”-二（1 - y2）y y = 0就是非线性的。2. 常微分方程的解及MATLAB旨令一阶常微分方程与高阶微分方程可以互化，已知一个n阶常微分方程（显式）y = f（t,y ,y , ,y ）（2）若令yy, yy, yyy（n J），可将上式化为n个一阶常微分方程组:5 = fi(t, yi, y2,yn)yi = f2(t, yi, y2,yn)(3)HI HIyn = fn(t, yi, y2,yn

4、)(3)式称为状态方程，y!, y2,yn (即y, y； y,y(n-1)称为状态变量，其中yi (即y)就是常微分方程（2）式的解。（3）式中右边的函数fi、f2、fn代表各个状态变量的一阶导 数的函数表达式，对于具体的方程它们有具体的形式，例如下列二阶非线性微分方程：y“_J(1_y2)y y =0若令yi二y, y2 = y，可将其改写成2个一阶微分方程组(状态方程)的形式：yi十W =(1 - y1)y2 - yi因此 fi *2； f2 -(i - y：)y2 - %。解析解只有少部分的线性常微分方程可以解析地求解(即可以算出精确的解表达式)，例如一阶常系数常微分方程 dy;dt

5、= y i可以通过直接积分解出，而多数微分方程尤其是非线性 方程则很难得到解析解。有解析解的方程虽然可以手算解出，但是MATLAE也提供了 dslove指令来求方程的解析 解，其使用格式：S = dsolve (方程1方程2；,初始条件1初始条件2，自变量) 方程用字符串表示，自变量缺省值为to i阶导数用D表示，2阶导数用D2表示，以此类推。 S用于返回方程解析解的表达式。如果是求解方程组，则S为一个结构体数组，它的每个域存放方程组每一个解的表达式。例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)%以最简形式显示 sans =-1/3*si n(x)*(-5+2*cos(x)% 方程的解(符号表达式)数值解对于没有解析解的方程主要依靠计算机进行数值求解(得到的是近似解)，例如方程y”.二(1 一 y2)y y =0就须通过计算机数值求解(结果是一系列解的数值而非表达式)。考虑n阶微分方程(2)式的数值求解，它等价于一阶常微分方程组(3)式。现将(3)式写成矩阵形式：Y(t)二 F(t,Y(t)其中yi(t)y、Y(t)=y2(t)*=y*(5) 为n个分量的列向量(Column vector)，也称状态向量y(t)创丿

7、f2(6)n个分量的列向量，其每个兀素分别为(3)式右边的函数表达式我们知道，微分方程要有唯一的确定解, 必须给定初值条件(to为初始时刻)：必须给定初值条件。因此方程(4)式要有确定的解r y(to)、 y(to)yi(to)Ay2(t。)丫0 =丫 (t0)=*=*yn (to)(7)Matlab提供了 ode45指令(ode是常微分方程的英文缩写)来求解方程(4)的数值解(近似解)。基本使用格式：tout, Yout =ode45 (odefun, tspan, Y0,options)其参数说明如下：odefun 般是用 M文件编写的函数，odefun代表函数名，由用户自己定义。函数返

8、回值为 式右边的F(t,Y)= (f1, f2，,fn) T。故odefun函数的返回值应是列向量，其最简单 的编写格式为：其中 tfunction F = odefun (t, Y) 【作用是计算并返回4式中的F(t, Y)】时间变量，为标量，代表计算进程中的某时刻点代表状态变量的列向量(即5式)返回值F(t, Y)，为列向量(见 6 式)ode45求解指令在计算时将会不断地在各个时间点调用odefun函数，并自动给输入参数t和Y赋值。tspan指定方程的求解区间t, tf，t0是初始时刻。Y0用户给定的初值条件，为n个分量的列向量，见(7)式。options可选项。一般情况下可缺省即可，若

9、用户有特殊要求则须使用odeset指令设置options选项，具体用法可使用help odeset命令查询，此处不做要求。tout列向量，输出求解过程中区间如tf上各个计算点的时刻，即tout = (to t t| tf ) , to, tf上计算点的数目是由Matlab自动生成。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时的零输入响应，就必须求解(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为字符型数组(字串)，其值为方程9的解表达式L=0.5;C=0.02;R=12.5;% 元件参数t=0:0.01:1;%定义区间0,1上的时间序列y=eval(S);% eval串演算函数，计算字串S (即9式的解)在t时刻的值plot(t,y)%绘制电压波形，即uc(t)的零输入响应B.数值解方程(9)虽然可以解析求解，但也可以使用则方程9改写成如下状态方程：ode45指令来近似计算。令 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)式的形式，因此状态向量第一步：通过 M文件创建 ode函数，函数名 circuit_2order_odefun。functionF=circuit_2order_odefu n( t,Y)global L C R %定义全局变量 L、C、R,以实现参数在 MATLAB的基本工作空 %间和函数的专用空间之间数据的传递F=Y(2);-1/L/C*Y(1)-R/L*Y(2);% 函数返回值 F,列向量，见(12 )式注意，函数创建完后，必

12、须存盘，存储文件名和函数名应一致！第二步：创建 M脚本文件circuit_2order.mL=0.5; C=0.02; R=12.5;% 元件参数Y0=1;0;%列向量，初值条件tspan=0,1;%定义求解区间0,1tout,Yout=ode45(circuit_2order_odefu n,tspa n,Y0);plot(tout,Yout(:,1);%绘电压波形，即uc(t)的零输入响应2. 二阶非线性电路一范德堡(Van de Pol)电路范德堡电路由一个线性电感、一个线性电容和一个非线性电阻构成，如图(a)所示。非线性电阻的伏安特性如图(b)所示，可很明显看出电阻是非线性的。范德堡电

13、路的特性可由一个二阶非线性微分方程描述(参见邱关源的电路第5版P463-P465):d2 iLdt2其中；-CT。上式就是著名的范德堡方程。该方程很难解析求解， 必须借助计算机数值求解。令力=畀2二di dt,将范德堡方程改写成一阶微分方程组（状态方程）的形式：= y2y2(1 -y：)y2 - y(14)假定参数；=0.1，求解区间上式也可以写成 式的形式。因此状态向量 Y和F分别为5(0)、2（0）J丿丫0(16)75匸(f 13、Y =|，F = | =2也丿lf2丿严(1 y ) y2 - y1 丿(15)0, 100，以及初值条件请仿照实验内容1,利用ode45指令求解范德堡方程，并

14、绘制电流Il的波形（即电流的零输入响应图）。四、实验任务1. 输入实验内容1中提供的程序上机练习，熟悉dslove指令以及ode45指令的使用方法,并尝试解答思考题1。2. 在任务1的基础上，独立编制实验内容2的计算程序，并尝试解答思考题3。五、思考题1. 在实验内容1中，若使用ode45指令解方程9式，如何得到电流的零输入响应波形？2. 在实验内容1中，试研究R=12 1, 3 1,，10时电压Uc（t）的响应波形，并把它们 绘制在一张图上。3. 从任务2的电流波形图（如下图）中，你能得到什么结论？如何绘制范德堡电路的uc Il函数图像（称为相图）？（提示：u _丄diL . 1 -3 _!

15、 y .1 y3 y ）UcIIIIy1 _ ys dt 3名 3附：实验任务2的参考图实验二思考题 1函数 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報国 CMend3任务二函数 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