系统建模与仿真 高阶电路分析

1、作业1：蒲丰投针试验模拟一实验原理： 根据频率的稳定性，当投针次数n很大时，测出针与平行线相交的次数m,则频率值m/n可作为P(A)的近似值代入上式，那么 利用上式可计算圆周率的近似值。二实验模拟过程 可以采用MATLAB软件进行模拟实验，即用MATLAB编写程序来进行“蒲丰投针实验”，并用GUI图形界面化直观显示。 Matlab的GUI编程代码如下： clcclose allclear allh0=figure('toolbar','none','position',350 180 570 530,'color',0.5 0.8

2、 0.9,'name','布丰投针试验');h_panel=uipanel('Units','normalized','position',0.05 0.72 0.35 0.25,'title','输入参数：','fontsize',11,'FontWeight','bold');h_text1=uicontrol(h_panel,'style','text','Units','n

3、ormalized','HorizontalAlignment','left','position',0.05 0.62 0.5 0.25,'string','两线间的宽度 a：','fontsize',10);h_edit1=uicontrol(h_panel,'style','edit','Units','normalized','HorizontalAlignment','left',

4、9;position',0.60 0.68 0.2 0.18,'fontsize',10);h_text2=uicontrol(h_panel,'style','text','Units','normalized','HorizontalAlignment','left','position',0.05 0.42 0.6 0.25,'string','针的长度 L:','fontsize',10);h_edit2

5、=uicontrol(h_panel,'style','edit','Units','normalized','HorizontalAlignment','left','position',0.60 0.48 0.2 0.18,'fontsize',10);h_text3=uicontrol(h_panel,'style','text','Units','normalized','Horizont

6、alAlignment','left','position',0.05 0.22 0.6 0.25,'string','投针次数 N：','fontsize',10);h_edit3=uicontrol(h_panel,'style','edit','Units','normalized','HorizontalAlignment','left','position',0.60 0.28 0.2

7、0.18,'fontsize',10);h_button1=uicontrol('Units','normalized','position',0.03 0.50 0.18 0.08,'string','开始投针','fontsize',10,'callback','function2');h_button3=uicontrol('Units','normalized','position',0.21 0

8、.50 0.18 0.08,'string','暂停','fontsize',10,'callback','function3');h_button2=uicontrol('Units','normalized','position',0.05 0.30 0.30 0.06,'string','统计','fontsize',10,'callback','function1');h_pane

9、l1=uipanel('Units','normalized','position',0.05 0.10 0.8 0.15,'title','数据分析：','fontsize',11,'FontWeight','bold');h_edit4=uicontrol(h_panel1,'style','edit','Units','normalized','HorizontalAlignment'

10、;,'right','position',0.35 0.05 0.6 0.35,'fontsize',10);h_text4=uicontrol(h_panel1,'style','text','Units','normalized','HorizontalAlignment','left','position',0.05 0.12 0.3 0.25,'string','估算的结果 ：','fon

11、tsize',10);h_edit5=uicontrol(h_panel1,'style','edit','Units','normalized','HorizontalAlignment','right','position',0.35 0.55 0.2 0.35,'fontsize',10);h_text5=uicontrol(h_panel1,'style','text','Units','norm

12、alized','HorizontalAlignment','left','position',0.05 0.55 0.3 0.25,'string','与平行线相交次数 y：','fontsize',10);h_text6=uicontrol('style','text','Units','normalized','HorizontalAlignment','left','positio

13、n',0.32 0.23 0.5 0.04,.'BackgroundColor',0.5 0.8 0.9,'FontWeight','bold','ForegroundColor',0.9 0 0,'string','应用近似公式: =(2*N*L)/(a*m)','fontsize',12);h_axes=axes('position',0.45 0.4 0.53 0.53,'YGrid','on','GridLine

14、Style','-');h_line=line('color',0 0.5 0.5,'linestyle','.','markersize',2,'erasemode','none'); %function1.ma=str2num(get(h_edit1,'string'); L=str2num(get(h_edit2,'string'); N=str2num(get(h_edit3,'string'); f=unifrnd(0

15、,pi,N,1); x=unifrnd(0,L,N,1); y=x<0.25*a*sin(f); m=sum(y); set(h_edit5,'string',num2str(m); format long; x1=2*N*L; x2=a*m; p=vpa(x1/x2); set(h_edit4,'string',num2str(double(p); %function2.mglobal k;while j<=100L=str2num(get(h_edit2,'string')*100;N=str2num(get(h_edit3,&#

16、39;string')*100;x=zeros(1,L);y=zeros(1,L);X=rand(1,N);Y=rand(1,N);angle=pi*rand(1,N);k=floor(N*rand(1,1);P=X(k);Q=Y(k);W=N;j=0;k=1;for i=1:L x(1)=P; y(1)=Q; set(h_line,'XData',x(i),'YData',y(i); x(i+1)=x(i)+0.001*cos(angle(k); y(i+1)=y(i)+0.001*sin(angle(k);end pause(0.005);j=j+1

17、;if (k=0) break;endend %function3.mglobal k;k=0; GUI界面如下：输入相关参数，并演示如下：逐渐增大N值N=10000时，PI=3.1625553447185326128021642944077N=100000时，PI=3.1467321186947354583196556632174N=1000000时，PI= 3.1459401642180764291367722762516三实验结论从上述数据分析可知，随着模拟次数的越来越多，PI的值逐渐稳定在值附近，即越来越趋近于，故蒲丰投针实验确实可以模拟出的值。作业2：高阶电路的系统建模与时域分析实验

18、原理： 高阶LC电路图如图所示： 对于一个N 阶线性电路，可以列出描述这个N阶电路的数学方程N 阶常微分方程。从此方程解得响应变量 ( )，这种分析计算动态电路的方法称为经典法。经典法是一种比较古老的方法，其优点是在分析低阶电路时，物理意义较清楚，容易掌握电路过度过程的规律。但用此法分析高阶电路时就要困难一些，计算工作量很大，这是因为 ( )的高阶导数的初始值常是未知的，需要从已知的其他初始条件导出。故采用状态方程分析方法，过程如下： 适当的选择电路中N个状态变量，可以将一个N阶微分方程变换为N个一阶微分方程，它表示N个状态变量间以及他们与激励间的关系，对于本题采用uc1,iL1,uc2,il

19、2,uc3,iL3作为中间状态变量，根据各回路电压电流关系得出如下方程： C3duc3dt=iL2(t)-iL3(t) (1)L3 diL3dt= uc3(t)+R32(i2(t)-i3(t)-R31iL3(t) (2)C2duc2dt=iL1(t)-iL2(t) (3)L2diL2dt= uc2t)+ R22(i1(t)-i2(t)-R31iL3(t)- iL2(t)R31- uc3(t)+R32(i2(t)-i3(t)-R31iL3(t) (4)C1duc1dt=i(t)-iL1(t) (5) L1diL1dt= uc1t)+ R12(i(t)-iL1(t)-R31iL3(t)- uc3(

20、t)+R32(i2(t)-i3(t)-R31iL3(t) R21Il2(t)- R11iL1(t)- uc2t)+ R22(i1(t)-i2(t)-R31iL3(t)- iL2(t)R31-uc3(t)+R32(i2(t)-i3(t)-R31iL3(t) (6) 化简为矩阵型式为： duc3dtdiL3dtduc2dtdiL2dtduc1dtdiL1dt=01L30-1L200 -1C3-R31+R32L30R32L200 0001L20-1L1 1C3R32L3-1C2-R22+R21+R320R22L1 000001L1 001C2R22L2-1C1-R22+R21+R11L1uc3iL3

21、uc2il2uc1iL1+00001/C1R21/L1i(t) 状态方程的求解： 状态方程的一般形式为 X=AX+B·S为了解此维常微分状态方程，使用 Matlab专门用于解微分方程的功能函数ode23函数。编写的Matlab程序源代码如下： % 清除内存、变量以及关闭Figureclc close allclear all% Global variable, 定义全局变量global iis;global r11;global r12;global l1;global c1;global r21;global r22;global l2;global c2;global r31;g

22、lobal r32;global l3;global c3;r31 = 5e-3;r32 = 0.5e-3;l3 = 10e-9;c3 = 30e-6;r21 = 5.5e-3;r22 = 0.8e-3;l2 = 1.3e-9;c2 = 1e-6;r11 = 10.8e-3;r12 = 5e-3;l1 = 103e-12;c1 = 30e-9;% time step for simulationt_step = 1e-10;t_end = 1e-06;tt = 0:t_step:t_end+t_step;% Create space to store temporary variablevtn

23、oise1 = zeros(1,floor(t_end/t_step)+1); % voltage of the branch of networkvcnoise1 = zeros(1,floor(t_end/t_step)+1); % voltage of the capacitor1icnoise1 = zeros(1,floor(t_end/t_step)+1); % current of the capacitor1ilnoise1 = zeros(1,floor(t_end/t_step)+1); % current of the inductance 1vcnoise2 = zer

24、os(1,floor(t_end/t_step)+1); % voltage of the capacitor2ilnoise2 = zeros(1,floor(t_end/t_step)+1); % current of the inductance 2vcnoise3 = zeros(1,floor(t_end/t_step)+1); % voltage of the capacitor3ilnoise3 = zeros(1,floor(t_end/t_step)+1); % current of the inductance 3% Initialize the indextotalind

25、ex = 1;peaktopeak = zeros(1, floor(t_end/t_step)+1);for iin=1:length(tt)-2 iis=0.5; tv,Yv=ode23('funsys2stage', 0 t_step, vcnoise1(iin); ilnoise1(iin);vcnoise2(iin); ilnoise2(iin);vcnoise3(iin); ilnoise3(iin); nYv=length(Yv(:,1); vcnoise1(iin+1)=Yv(nYv,1); ilnoise1(iin+1)=Yv(nYv,2); icnoise1

26、(iin+1) = iis-ilnoise1(iin+1); vtnoise1(iin+1) = vcnoise1(iin+1)+r12*icnoise1(iin+1);end% Update the variablevtnoise1 = vtnoise1(1,2:length(vtnoise1);peak = max(vtnoise1); valley = min(vtnoise1);peaktopeak(1,totalindex) = peak-valley;totalindex = totalindex+1;%查找u(t)的第一零点位置，从而确定最差情况下方波的频率m = length(vtnoise1); x1=vtnoise1(1:m-1); x2=vtnoise1(2:m); indz = find(vtnois