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

clear, format compact
L=0.5; C=0.02; R=12.5;		% 輸入元件參數
uc0=1; iL0=0;
for R=1:10
   a=R/2/L;; w0=sqrt(1/(L*C)); % 輸入給定參數
   a1=a-sqrt(a^2-w0^2);		% 方程的兩個根
   a2=a+sqrt(a^2-w0^2);
   dt=0.01;t=0:dt:1;				% 設定時間數組
   % 方法1，用公式
   uc1=(a2*uc0+iL0/C)/(a2-a1)*exp(-a1*t);		% uc的第一個分量
   uc2=-(a1*uc0+iL0/C)/(a2-a1)*exp(-a2*t);		% uc的第二個分量
   iL1=-a1*C*(a2*uc0+iL0/C)/(a2-a1)*exp(-a1*t);		% iL的第一個分量
   iL2=a2*C*(a1*uc0+iL0/C)/(a2-a1)*exp(-a2*t);		% iL的第二個分量
   uc=uc1+uc2; iL=iL1+iL2;							% 把兩個分量相加
   % 分別畫出兩種數據曲線
   figure(1),plot(t,uc),hold on
   figure(2),plot(t,iL),hold on
end
   figure(1),grid,figure(2),grid
% 方法2，用拉普拉斯變換及留數法
for R=1:9
   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=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t);	% 求時域函數
   iLn=C*diff(ucn)/dt;			% 對ucn求導得到電流iLn
   figure(3),plot(t,ucn),grid,hold on		% 繪曲線
   figure(4),plot(t(2:end),iLn),grid,hold on
end 
figure(3),grid,figure(4),grid