%《MATLAB在電子信息課程中的應用》第五章例5.6程序q506
% 直流電路的暫態計算：二階過阻尼電路
% 電子工業出版社出版  陳懷琛 吳大正 高西全合著  2001年10月

clear, format compact
L=0.5; R=12.5; C=0.02;		% 輸入元件參數
uc0=1; iL0=0;
alpha=R/2/L; wn=sqrt(1/(L*C)); % 輸入給定參數,這里ρ即α。
p1=-alpha+sqrt(alpha^2-wn^2)		 	  % 方程的兩個根
p2=-alpha-sqrt(alpha^2-wn^2)
dt=0.01;t=0:dt:1;			  % 設定時間數組
% 方法1,用公式
uc1=(p2*uc0-iL0/C)/(p2-p1)*exp(p1*t);		% uc的第一個分量
uc2=-(p1*uc0-iL0/C)/(p2-p1)*exp(p2*t);		% uc的第二個分量
iL1=p1*C*(p2*uc0-iL0/C)/(p2-p1)*exp(p1*t);		% iL的第一個分量
iL2=-p2*C*(p1*uc0-iL0/C)/(p2-p1)*exp(p2*t);		% iL的第二個分量
uc=uc1+uc2; iL=iL1+iL2;									% 把兩個分量相加
% 分別畫出兩種數據曲線
subplot(2,1,1),plot(t,uc),grid
subplot(2,1,2),plot(t,iL),grid
disp('方法2,用拉普拉斯變換及留數法')
pause
num=[uc0,R/L*uc0+iL0/C];		% uc(s)的分子系數多項式
den=[1,R/L,1/L/C];				% uc(s)的分母系數多項式
[r,p,k]=residue(num,den);		% 求極點留數
% 求時域函數ucn, 對ucn求導得到電流iLn
ucn=r(1)*exp(p(1)*t)+r(2)*exp(p(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),grid