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