MATLAB在电路中的应用

MATLAB应用（三）——Matlab在电路中的应用2

MATLAB中的变量与常量都是矩阵(标量可看做1×1阶的矩阵，向量可看做n×1或1×n阶的矩阵)，其元素可以是复数和任意形式的表达式，它具有元素群运算能力。MATLAB的这些优于其他语言的特色，有利于分析计算电路的各种问题，并且使编程更简便，运算效率更高。3

学习目的：

通过介绍计算电路问题的编程方法和技巧，逐步熟悉MATLAB语言的使用。

例题的解法本身，不一定最佳。

求解电路的专用软件：Spice、PSpice等软件4内容:

电阻电路的求解(例1-3)动态电路的求解(例4-7)例题分析过程：

例题说明求解过程：建模Matlab程序说明Matlab程序运行、结果演示5电阻电路的求解

如us=10V,

求i3，u4，u7；(2)

如已知u4=6V，

求us,i3,u7。图1例1的电阻电路

[例1]

如图1所示的电路，己知：R1=2Ω，R2=4Ω，R3=12Ω,

R4=4Ω,R5=12Ω,R6=4Ω,R7=2Ω。6

对图示电路，用网孔电流法列写网孔电流方程如下：建模解：7写成矩阵形式为：也可直接列写数字方程为：R1=2Ω，R2=4Ω，R3=12Ω，R4=4ΩR5=12Ω，R6=4Ω，R7=2Ω8

矩阵方程简写为：

令us=10V，求解矩阵方程得到ia、ib、ic。再由i3=ia-ib

，u4=R4ib

，u7=R7ic即可得到问题(1)的解9

根据电路的线性性质，可令i3=k1us，u4=k2us,u7=k3us,由问题(1)的解求得比例系数，进一步使问题(2)得到解答。具体根据问题(1)的结果可列出以下的表达式：因此，通过下列表达式即可求得问题(2)的解：10Matlab程序(Ex01.m)clear,closeall,

formatcompactR1=2;R2=4;R3=12;R4=4;R5=12;R6=4;R7=2;

%为给定元件赋值display('解问题(1)')

%解问题(1)a11=R1+R2+R3;

a12=-R3;

a13=0;

%将系数矩阵各元素赋值a21=-R3;

a22=R3+R4+R5;

a23=-R5;a31=0;

a32=-R5;

a33=R5+R6+R7;b1=1;b2=0;b3=0;us=input(‘

给定us=’),

%输入解(1)的已知条件A=[a11,a12,a13;a21,a22,a23;a31,a32,a33];

%列出系数矩阵AB=[b1;0;0];I=A\B*us;%

I=[ia;ib;ic]ia=I(1);

ib=I(2);

ic=I(3);i3=ia-ib,u4=R4*ib,u7=R7*ic%解出所需交量display('解问题(2)')

%利用电路的线性性质及问题(1)的解求解问题(2)u42=input('给定u42=');k1=i3/us;

k2=u4/us;

k3=u7/us;

%由问题(1)得出待求量与us的比例系数us2=u42/k2,

i32=k1/k2*u42,

u72=k3/k2*u42

%按比例方法求出所需交量11程序运行结果解问题(1)给定us=10i3=0.3704，u4=2.2222，u7=0.7407解问题(2)给定u42=6us2=27.0000，i32=1.0000，u72=2运行结果：电路的解：i3=0.3704A,

u4=2.2222V,

u7=0.7404Vus=27V,

i3=1A,

u7=2VEx01.m12补充说明：

实际中，如果熟悉列方程的方法,那么在编写MATLAB程序时可直接写出A和B为：从而可省去给元件和矩阵各元素赋值等语句。13[例2]

对如图2所示的电路，已知R1=R2=R3=4Ω,R4=2Ω,控制常数K1=0.5,k2=4,is=2A,求i1和i2。图2例2的电路14

对图示电路，用节点电压法列写方程得：建模解：uaub15

根据图示电路，控制变量i1、i2与节点电压ua、ub的关系为：

整理以上两式，将i1、i2也作为未知量，和前面的节点电压共同组成方程

,并写成矩阵形式有：

令is=2A，求解上式即可得到i1和i2

。uaub16Matlab程序(Ex02.m)clear,

formatcompactR1=4;R2=4;R3=4;R4=2;%设置元件参数is=2;k1=0.5;k2=4;%按A*X=B*is列写电路的矩阵方程，其中X=[ua;ub;i1;i2]。a11=1/R1+1/R2;a12=-1/R2;a13=0;a14=-k1;%设置系数Aa21=-1/R2;a22=1/R2+1/R3+1/R4;a23=-k2/R3;a24=k1;a31=1/R2;a32=-1/R2;a33=-1;a34=0;a41=0;a42=1/R4;a43=0;a44=-1;A=[a11,a12,a13,a14;a21,a22,a23,a24;a31,a32,a33,a34;a41,a42,a43,a44];B=[1;0;0;0];%设置系数BX=A\B*is;i1=X(3),i2=X(4)%显示要求的分量

17Ex02.m程序运行结果(电路的解)i1=1，i2=118

[例3]

对如图3所示的电路，已知R1=4Ω,R2=2Ω,R3=4Ω，R4=8Ω；is1=2A，is2=0.5A。图3例3的电路(1)负载RL为何值时能获得最

大功率?(2)

研究RL在0~10Ω范围内变

化时,其吸收功率的情况。19解：

用戴维南等效电路来求解。对图3(a)电路，断开ao,并在ao端接入外电流源ia,

如图3(b)所示。以o为参考点列节点方程得：建模图3例3的电路20前面的方程写成矩阵形式为：其中：戴维南等效电路如图3(c)所示，其方程为：图3例3的等效电路21方法Ⅰ：

令ia=0,

is1=2A,

is2=0.5A,

由矩阵方程求得u11，u21，ua1。因ia=0,由戴维南等效电路方程得：uoc=ua1。

再令is1=is2=0，ia=1A,仍由矩阵方程可求得另一组u12，u22，ua2。由于内部电源is1=is2=0，故uoc=0。从而由戴维南等效电路方程有：22

于是，原电路戴维南等效电路如图3(d)所示，负载RL获得最大功率时有:图3例3的等效电路

至于问题(2),由图3(d)可得RL吸收功率为：再令RL=lΩ，2Ω，3Ω，…，1OΩ，即可由上式分别求得PL,

并画图。23

可设ia为一个序列(如ia=0.1，0.2，…，2)，计算相应的ua序列，再用线性拟合，得出如下的直线方程：方法Ⅱ：

从而求得：24Matlab程序Ⅰ(Ex03-1.m)clear,

formatcompactR1=4;

R2=2;

R3=4;

R4=8;

%设置元件参数is1=2;is2=0.5;%按A*X=B*is列写此电路的矩阵方程,其中X=[u1;u2;ua];is=[is1;is2;ia]a11=1/R1+1/R4;a12=-1/R1;a13=-1/R4;%设置系数矩阵a21=-1/R1;a22=1/R1+1/R2+1/R3;a23=-1/R3;a31=-1/R4;a32=-1/R3;a33=1/R3+1/R4;A=[a11,a12,a13;a21,a22,a23;a31,a32,a33];B=[1,1,0;0,0,0;0,-1,1];%设置系数矩阵B%方法Ⅰ:令ia=0,求uoc=x1(3);再令is1=is2=0,设ia=1,求Req=ua/ia=x2(3).Xl=A\B*[is1;is2;0];uoc=X1(3)X2=A\B*[0;0;1];Req=X2(3)RL=Req;P=uoc^2*RL/(Req+RL)^2%求最大负载功率%也可设RL为一数组,求出的负载功率也为一数组,画出曲线找极大值RL=0:10,p=(RL*uoc./(Req+RL)).*uoc./(Req+RL),%设RL序列,求其功率figure(1),plot(RL,p),grid %画出功耗随RL变化的曲线xlabel(‘RL'),ylabel(‘p')25Matlab程序Ⅱ(Ex03-2.m)clear,

formatcompactR1=4;

R2=2;

R3=4;

R4=8;

%设置元件参数is1=2;is2=0.5;%按A*X=B*is列写此电路的矩阵方程,其中X=[u1;u2;ua];is=[is1;is2;ia]a11=1/R1+1/R4;a12=-1/R1;a13=-1/R4;%设置系数矩阵a21=-1/R1;a22=1/R1+1/R2+1/R3;a23=-1/R3;a31=-1/R4;a32=-1/R3;a33=1/R3+1/R4;A=[a11,a12,a13;a21,a22,a23;a31,a32,a33];B=[1,1,0;0,0,0;0,-1,1];%设置系数矩阵B%方法Ⅱ:设一个ia序列,计算一个ua序列,用线性拟合求出其等效开路电压和等效内阻fork=1:21ia(k)=(k-1)*0.1;X=A\B*[is1;is2;ia(k)];%定义X=[u1;u2;ua]u(k)=X(3);endfigure(2),plot(ia,u,’x’),grid%线性拟合xlabel('ia'),ylabel('ua')c=polyfit(ia,u,1);%ua=c(2)*ia+c(1),用拟合函数求c(1),c(2)uoc=c(1),Req=c(2)26程序运行结果uoc=5V,

Req=5Ω,

Pmax=1.25W.(a)功率随负载的变化曲线Ex03_1.M27程序运行结果(b)电路对负载的输出特性Ex03_2.M28动态电路的求解[例4]

一阶动态电路如图4所示，己知:Rl=3Ω,R2=2Ω,R3=6Ω，C=1F；us=18V，is=3A，在t<0时，开关S位于“1”，电路已处于稳定状态。图4动态电路(1)t=0时，开关S闭合到“2”，

求uc，iR2(t)，并画出波形；(2)若经10秒，开关S又复位到“1”，

求uc(t)，iR2(t)，并画出波形。29对该一阶动态电路可用通用的解决方案[式(2.33)](也称三要素法)求解。建模解：(2.33)

首先求初始值uc(O+)和iR2(O+)。

为此，先求uc(O-)，在t=0-时，开关位于“1”，电路已达到稳定。电容可看做开路，不难求得uc(O-)=-12V。30

根据换路定则(电容电压不变)，得电容初始电压uc(O+)=uc(O-)=-12V。

在t=0时，开关己闭合到“2”，可求得非独立初始值iR2(O+)为：31

其次求稳定值。

达到稳态时电容可看做开路，于是可得：

时间常数为：

因此，解为：(2.33)32

经10秒后，开关又闭合到“1”，将t=10代入前面的电压表达式可得电

容电压的初始值为：

由图可见这时

并保持不变。

达到稳定时，这时时间常数为：

33

利用通用公式，得到uc(t)、iR2(t)为：(2.33)34Matlab程序(Ex04.m)clear,closeall,

formatcompactR1=3;

us=18;

is=3;

R2=12;

R3=6;

C=1;

%给出原始数据%解问题(1)uc0=-12;ir20=uc0/R2;ir30=uc0/R3;%算出初值ir20及ucOic0=is-ir20-ir30;ir2f=is*R3/(R2+R3);%算出终值ir2f及ucfir3f=is*R2/(R2+R3);ucf=ir2f*R2;icf=0;%注意时间数组的设置,在t=O及10附近设两个点,见图4(a)t=[[-2:0]-eps,0:9,10-eps,10+eps,11:20];figure(1),plot(t),grid%从图(a)中可看出时间与时间数纽下标的关系,t=10+eps对应下标15uc(1:3)=-12;ir2(1:3)=3;%t<0时的值T=R2*R3/(R2+R3)*C;%求充电时常数uc(4:14)=ucf+(ucO-ucf)*exp(-t(4:14)/T);%ir2(4:14)=ir2f+(ir20-ir2f)*exp(-t(4:14)/T);%用三要素法求输出35%解问题(2)uc(15)=uc(14);ir2(15)=is;%

求t=10+eps时的各初值ucf2=-12;

ir2f=is;%求uc和ir2在新区间终值ucf2和ir2fT2=R1*R3/(R1+R3)*C;

%再用三要素法求输出uc(15:25)=ucf2+(uc(15)-ucf2)*exp(-(t(15:25)-t(15))/T2);ir2(15:25)=is;figure(2)subplot(2,1,1);h1=Plot(t,uc);

%绘uc图grid,

set(h1,‘linewidth’,2)

%加大线宽subplot(2,1,2),

h2=plot(t,ir2);

%绘ir2图grid,

set(h2,'linewidth',2)

36程序运行结果(a)时间与其数组下标的关系(b)uc及iR2的暂态波形

Ex04.m37[例5]

如图5所示的一阶电路,已知R=2Ω，C=0.5F,电容初始电压uc(O+)=4V,激励的正弦电压us(t)=umcosωt,

其中um=10V，ω=22rad/s。图5正弦激励一阶电路

当t=0时，开关S闭合，求电容电压的全响应，区分其暂态响应与稳态响应,并画出波形。38

电路中电容电压的微分方程为：建模解：其时间常数为：39(2.33)

用三要素法，其解为：式中uc(O+)为电容初始电压，ucp(t)为微分方程的特解。

当正弦激励时，设ucp(t)=ucmcos(ωt+φ),其中：40最后得电容电压的全响应为：其暂态响应(固有响应)为：稳态响应(强迫响应)为:41Matlab程序(Ex05.m)clear,closeall,

formatcompactR=2;

C=0.5;

T=R*C;

ucO=4;

%输入元件参数um=10;

w=2;

Zc=1/j/w/C;%输入给定参数t=0:0.1:10;

%设定时间数组us=um*cos(w*t);

%设定激励信号ucst=us*Zc/(R+Zc);

%稳态分量计算ucpO=ucst(1);

%稳态分量的初始值uctr=(ucO-ucpO)*exp(-t/T);

%暂态分量uc=uctr+ucst;

%总的uc为两项之和%把三种数据画在一张图上plot(t,

uc,

‘-’,t,uctr,

':',

t,

ucst,

‘-.'),

gridLegend(‘uc‘,'uctr‘,‘ucst’)

%用图例标注42程序运行结果电容上的电压波形Ex05.m43

[例6]

考察二阶过阻尼电路的固有响应(零输入响应),图6为典型的二阶动态电路，其固有响应有过阻尼、临界阻尼和欠阻尼三种情况。此例讨论过阻尼情况。图6过阻尼二阶电路

己知L=0.5H,C=0.02F，R=125Ω，初始值uc(O)=1v，iL(O)=0，求t≥0时的uc(t)和iL(t)的固有响应,并画出波形。44按图不难列出关于uc的微分方程为：建模解：方法Ⅰ：令衰减常数

，谐振角频率

，则得二阶微分方程为：45即α>ω0，表现为过阻尼，其解为：式中：在此，过阻尼解其初始值为：46

对微分方程作拉氏变换，考虑到初始条件，可得：方法Ⅱ：整理可得：对上式求拉氏反变换即可得到时域的表达式，将等式右端的多项式分解为部分分式，得：47其中num和den分别为分子、分母多项式系数组成的数组。进而写出：

s1，s2，r1和r2可以用代数方法求出,在MATLAB中有residue函数,专门用来求多项式分式的极点和留数，其格式为：这样就无需求出其显式，使得程序特别简明。上式中，sl和s2是多项式分式的极点，r1和r2是它们对应的留数。从而有：48Matlab程序(Ex06.m)clear,closeall,

formatcompactL=0.5;

R=12.5;

C=0.02;

%输入元件参数ucO=1;

iLO=0;alpha=R/2/L;

w0=sqrt(1/(L*C));

%输入给定参数s1=-alpha+sqrt(alpha^2-w0^2)%方程的两个根s2=-alpha-sqrt(alpha^2-w0^2)dt=0.01;t=0:dt:1;%设定时间数组%方法Ⅰ,用公式uc1=(s2*ucO-iLO/C)/(s2-sl)*exp(s1*t);%uc的第一个分量uc2=-(s1*ucO-iLO/C)/(s2-s1)*exp(s2*t);%uc的第二个分量iL1=s1*C*(s2*ucO-iLO/C)/(s2-sl)*exp(s1*t);%iL的第一个分量iL2=-s2*C*(s1*ucO-iLO/C)/(s2-s1)*exp(s2*t);%iL的第二个分量uc=uc1+uc2;iL=iL1+iL2;%把两个分量相加%分别画出两种数据曲线subplot(2,1,1),plot(t,uc),gridsubplot(2,1,2),plot(t,iL),grid49%方法Ⅱ,用拉普拉斯变换及留数法num=[uc0,R/L*ucO+iLO/C];%uc(s)的分子系数多项式den=[1,R/L,1/L/C];%uc(s)的分母系数多项式[r,s,k]=residue(num,den);

%求极点留数%求时域函数ucn,对ucn求导得到电流iLnucn=r(1)*exp(s(1)*t)+r(2)*exp(s(2)*t);iLn=C*diff(ucn)/dt;%%绘曲线,注意求导后数据长度减少一个figure(2),subplot(2,1,1),plot(t,ucn),grid%绘曲线subplot(2,1,2),plot(t(1:end-1),iLn),grid50程序运行结果电压uc和电流iL的波形Ex06_1.mEx06_2.mEx06.m51[例7]

考察二阶欠阻尼电路的固有响应(零输入响应)，电路同例6。如L=0.5H，C=0.02F。初始值uc(O)=lv,iL=0,试研究R分别为1Ω,2Ω,3Ω,…,1OΩ时,uc(t)和iL(t)的固有响应,并画出波形图。52电路的微分方程同例6,

为：建模解：其中

,谐振角频率

,且有。53在此ω0=10，当R=1Ω,2Ω,3Ω,…,1OΩ时，α=1,2,3,…

10，显然α=ω0=10为临界阻尼,其余为欠阻尼(衰减振荡)情况,这时方程的解为：式中同样可用拉氏变换及留数法求解，具体见程序。54Matlab程序(Ex07_1.m)clear,closeall,

formatcompactL=0.5;

C=0.02;

%输入元件参数ucO=1;

iLO=0;forR=1:10

alpha=R/2/L;

w0=sqrt(1/(L*C));

%输入给定参数

s1=-alpha-sqrt(alpha^2-w0^2);

%方程的两个根s2=-alpha+sqrt(alpha^2-w0^2);

dt=0.01;

t=0:dt:1;);

%设定时间数组

%方法1,用公式wd=sqrt(w0^2-alpha^2);

A=sqrt((wd*ucO)^2+(iLO/C+al