back
New '02/1/16 = 回路応答について別ファイルを作り始めた。
Update '01/6/29 = 学生用に用意されているソフトについて説明した。 Update '01/4/23 = Linkを設けた。(望月の他のページにあったものをここに再配置した） Update '00/7/9 = 一日体験入学に備えて言い回しを変更した。 Update '99/12/10 = my_func 周りの言い回しを訂正した。 Update '99/9/14 = 一部言い回しを訂正した。 Update '99/4/28 = 実行例を加え，一部言い回しを訂正した。 Update '98/12/11 = 微分方程式の項を訂正した(配列の使い方を変えた)。 Update '98/12/11 = 微分方程式の項を書き加えた。

MATLAB つかいませんか (File-2)
--- MATLAB を使うための ABC ---

応用２ （回路応答を MATLAB で調べる）
目次
  • 回路解析-1 −−− ２次の系の時間応答（微分方程式 ）
    • 組込み関数（ODE23など）を使って
    • 組込み関数を使わずに解く
    • Vr - Vout 平面上に解軌跡を表示する
  • 回路解析-2 −−− 周波数特性
    • 自分で多項式を計算
    • 組込み関数（polyval）を使って

回路解析-1 --- ２次の系の時間応答

回路の応答を調べることは，微分方程式を解くことと同じです。 ここでは，簡単な RLC 回路を例に取り，微分方程式を解きながら回路応答を探ります。
実行例 ( 他にも良いやり方はいろいろあります )
２次の系の時間応答 --- 微分方程式を解く
ここでは，LRC 直列回路に電源を加え， Cの両端から出力を取り出すことにします。 
回路図：
＋−−［Ｌ］−−−［Ｒ］−−＋ Vｏｕｔ
｜　　　　　　　　　　　　　｜
Ｖｉｎ　　　　　　　　　　［Ｃ］
｜　　　　　　　　　　　　　｜
---　　　　　　　　　　　　---
///　　　　　　　　　　　　///
方程式は，次式で与えられます。s は 時間微分を，1/s は積分を表わします。
    Vin = ( sL + R + 1/sC) i
    Vr = i R
    Vout = i/sC
なお，それぞれの定数は次のとおりとしておきます：
L=1mH, C=1 μF （これより， f0 は，約 5 kHz）， R=10 Ω 〜 100 Ω （これより，Q は大よそ 0.3〜3 です）

また，Vin は，正弦波 又は ステップ関数とします。

[解法-1 : 微分方程式として解く]
これを解くにあたって次の置き換えを考えます：
「2次の微分方程式は，2つの1次の微分方程式を連立させたものである」
s は微分を表わしていますので，上の式は次のように書き換えることができます。
    Vr ' = (- Vr - Vout + Vin) * R / L
    Vout' = Vr/(R * C)
まず，この関係式を関数にします。 下記プログラムのうち，xdot(1)は Vr ' を，xdot(2) は Vout' を表わします。 また，引数のうち，x(1) は Vr を， x(2) は Vout を表わします。 一般に，微分方程式は時間の関数でもあることから， 引数の中には時間 t も含まれます。 
 ---セーブするのは以下の6行 !! ファイル名は'my_func.m'--- 
function xdot=my_func(t,x)
 xdot=zeros(2,1); % 領域確保
 L=0.001; C=1.e-6; R=20;
 % Vin は T=1ms で立ち上がるステップ関数とする。
 xdot(1)=( -x(1) - x(2) + 0.5*(1+sign(t-0.001)) )*R/L;
 xdot(2)=x(1)/(R*C);
% end of function

% 例えば Vin が 300 Hz の正弦波ならば，xdot(1) の行を次行に変更。
xdot(1)=( -x(1) - x(2) + sin(2*pi*300*t) )*R/L;

MATLAB の組み込み関数 ode23 を使うと， 簡単に微分方程式を解けます。 早速"help ode23"と打ってマニュアルを確認しよう。 他にも，ode45のように，目的と使い方は同じだが， 内部処理（アルゴリズム）が異なる関数が用意されています。 
% 回路解析-1
% 時間領域: 0 --- 5ms (秒)
% xx(:,1) ... Vr   ;   初期条件は Vr = 0
% xx(:,2) ... Vout ;   初期条件は Vout = -1

[t,xx]=ode23('my_func',[0 0.005], [0 0]);
plot(t,xx(:,1),'b',t,xx(:,2),'r')

% end
[解法-2 : 組み込み関数を使わずに解く]
基本的には 解法-1 と同じアルゴリズムだが， 組み込み関数を使わず非常にシンプルに取り扱います。 （従って誤差は保証されない）
なお，解法-1 の my_func を流用します。 
% 回路解析-2

% 計算結果を入れる領域の確保
tt=zeros(10001,1);
vr=tt;
vout=tt;

% initial condition
dt=0.005/(length(tt)-1); % 時間: 0〜5 ms
tt(1)=0;
vr(1)=0; % 初期電流
vout(1)=0; % 初期出力電圧

% calculation
for k=1:length(tt)-1
 xdot=my_func(tt(k),[vr(k) vout(k)]);
 tt(k+1)  =  tt(k)+dt;
 vr(k+1) = vr(k)+dt*xdot(1);
 vout(k+1)=vout(k)+dt*xdot(2);
end

% display
plot(tt,vr,'b',tt,vout,'r')
% end

[解法-3 : Vr - Vout 平面上に解軌跡を表示する]
ここでも，もとの方程式を取り扱います。 解法-1 の my_func をそのまま流用します。 
% 解法-3 ： Vr - Vout 平面上での解軌跡
%	縦軸は x1( Vr ), 横軸は x2(Vout) 
[vvx,iiy]=meshgrid(-2:0.4:2, -0.8:0.2:0.8);
siz=size(vvx)
dx0=zeros(siz(1), siz(2), 2); dx1=dx0;
for m=1:siz(1)      % 各座標での解の動き方を表示します
   for n=1:siz(2)
      dx0(m,n,:)=my_func(0,[iiy(m,n) vvx(m,n)]);
      dx1(m,n,:)=my_func(1,[iiy(m,n) vvx(m,n)]);
   end
end
figure(1); quiver(vvx,iiy,dx0(:,:,2),dx0(:,:,1),'.')
hold on;   quiver(vvx,iiy,dx1(:,:,2),dx1(:,:,1),'g.')

% 回路解析3-1 : 時間領域 0 - 5ms，初期電流 0，初期電圧 -1
[t,xx]=ode23('my_func',[0 0.005], [0 -1]);
plot(xx(:,2), xx(:,1),'b')

% 回路解析3-2
tt=zeros(10001,1);    % 計算結果を入れる領域の確保
vr=tt; vout=tt;

dt=0.005/(length(tt)-1) % 時間領域: 0 〜 5 ms
tt(1)=0;
vr(1)=0;    % 初期電流
vout(1)=-1; % 初期出力電圧

for k=1:length(tt)-1
   xdot=my_func(tt(k),[vr(k) vout(k)]);
   tt(k+1) = tt(k)+dt;
   vr(k+1)  =  vr(k)+dt*xdot(1);
   vout(k+1)=vout(k)+dt*xdot(2);
end
plot(vout, vr, 'r'); hold off
figure(2); plot(t,xx(:,2),'b', tt,vout,'r')
% end


回路解析-2 --- 周波数特性

先ほどと同じ回路について，周波数応答を調べます。
実行例 ( 他にも良いやり方はいろいろあります )
周波数特性を求める
先ほどと同じ回路です。 LRC 直列回路に電源を加え， Cの両端から出力を取り出します。 
回路図：
＋−−［Ｌ］−−−［Ｒ］−−＋ Vｏｕｔ
｜　　　　　　　　　　　　　｜
Ｖｉｎ　　　　　　　　　　［Ｃ］
｜　　　　　　　　　　　　　｜
---　　　　　　　　　　　　---
///　　　　　　　　　　　　///
利得は，
    Vout／Vin = 1 ／ ( 1 + sCR + s2LC)
なお，それぞれの定数は次のとおりとしておきます：
L=1mH, C=1 μF （これより， f0 は，約 5 kHz）， R=10 Ω 〜 100 Ω （これより，Q は大よそ 0.3〜3 です）

[解法-1]
MATLAB は複素数を取り扱えますから，簡単に利得を計算できます。 とにかくそのまま書いていくだけです。 ただし，マトリクス同士の計算の乗除算は，".*", "./" です。 
% 回路解析-1
% frequency responce of an L-R-C circuit
%    各部の電圧の周波数依存性を求めるプログラム」

L=0.001; C=1.e-6; R=20;
f=logspace(1,5,100); s=i*2*pi*f;

gain=1./(1 + s.*C.*R + s.*s.*L.*C);

figure(1); loglog(f,abs(gain),'r')
figure(2); semilogx(f,angle(gain)*180/pi,'r')

% end

[解法-2]
polyval は多項式を取り扱う関数です。 計算時に，複素数も取り扱えますから，利得計算に使えます。 
% 回路解析-2
% frequency responce of an L-R-C circuit
%    各部の電圧の周波数依存性を求めるプログラム」

L=0.001; C=1.e-6; R=20;
f=logspace(1,5,100); s=i*2*pi*f;

gain=polyval([1],s) ./ polyval([L*C C*R 1],s);
%  分子 1 ， 分母 (s^2)LC + sCR + 1

figure(1); loglog(f,abs(gain),'r')
figure(2); semilogx(f,angle(gain)*180/pi,'r')

% end

//