ステップ応答の解析解

2.  数値計算

  （オイラー法と Runge-Kutta 法の精度比較）

3.  周波数応答とボード線図

•  ２次振動系と共振

1.  ステップ応答の解析解

2.  数値計算

3.  周波数応答とボード線図

•  レポート課題

2 次振動系の状態空間表現

€

d

²

x t ( )

dt

²

+ 2 ζω

_n²

^{dx t} ( )

dt + ω

_n²

^{x t} ( ) ^{= ω}

ⁿ²

^{u t} ( )

€

dx t ( )

dt = y t ( ) ^とおくと

€

d

²

x t ( )

dt

²

= dy t ( )

dt = − 2 ζω

_n²

^{y t} ( ) ^{− ω}

ⁿ²

^{x t} ( ) ^{+ ω}

ⁿ²

^{u t} ( )

よって、以下の状態空間表現が得られる

€

d dt

x y



 



  = 0 1

− ω

_n²

− 2 ζω

_n²



 



  x y



 



  + 0 ω

_n²



 



  u t ( )

€

x = [ 1 0 ] ^x

y



 



 

配布プログラム

•  ２次の Runge-Kutta 法 RK2.m

•  ベクトル場 SOE.m

•  メイン部分  SOE_main.m  

ベクトル場の関数 (SOE.m)

function xy = SOE(t,y,para)

   xy = [y(2);-para(2).^2.*y(1)-2.*para(1).*para(2).*y(2)+para(2).^2.*para(3)];

end

€

d dt

x y



 



  =

縦ベクトル

ステップ入力の大きさ

€

y

− ω

_n²

^x − 2 ζω

_n²

^y + ω

_n²

^{u t} ( )



  

 

時間 状態ベクトル パラメータ 速度ベクトル

function main

xi = 0.2; %減衰係数 wn = 2; %固有振動数

In = 1;%ステップ入力の大きさ dt = 0.1;%時間刻み

T = 15;%最終時間

xinit =[0;0];%初期値セット

para = [xi,wn,In];%パラメータセット

%@SOEで関数のポインタの受け渡し

[t1,xy] = RK2(@SOE,[0,T],xinit,dt,para); %2次Runge-Kutta x = xy(1,:); xmin = min(x)-0.1; xmax = max(x)+0.1;

y = xy(2,:); ymin = min(y)-0.1; ymax = max(y)+0.1;

subplot(2,1,1); %2つのグラフを縦に並べて表示 plot(t1,x);%ステップ応答の時間発展をプロット axis([0 T xmin xmax]);

xlabel('t'); ylabel('x');

subplot(2,1,2); %2つのグラフを縦に並べて表示 plot(x,y);%状態空間プロット

axis([xmin xmax ymin ymax]);

xlabel('x'); ylabel('dx/dt');

end

メイン部分

SOE_main.m

減衰係数を変えて、

ステップ応答の変化を 観察

MATLAB 演習（線形システム）

•  微分方程式の数値解法（復習）

•  １次遅れ系と時定数

1.  ステップ応答の解析解  2.  数値計算

  （オイラー法と Runge-Kutta 法の精度比較）

3.  周波数応答とボード線図

•  ２次振動系と共振

1.  ステップ応答の解析解 2.  数値計算

3.  周波数応答とボード線図

•  レポート課題

配布プログラム

•  ２次の Runge-Kutta 法 RK2.m

•  ベクトル場 SOE_FR.m

•  メイン部分  SOE_FR_main.m  

ベクトル場の関数 (SOE.m)

€

d dt

x y



  

  =

縦ベクトル

正弦波の振幅

€

y

− ω

_n²

^x − 2 ζω

_n²

^y + ω

_n²

^{A sin} ( ) ω ^t



 



 

時間 状態ベクトル パラメータ 速度ベクトル

function xy = SOE(t,y,para)

xy = [y(2);… %改行する場合、…

-para(2).^2.*y(1)-2.*para(1).*para(2).*y(2)+para(2).^2.*para(3).*sin(para(4).*t)];

end

メイン部分 SOE_FR_mai

n.m

減衰係数を変えて、

応答の振幅変化を観 察

xi = 0.5; %減衰係数 wn = 1; %固有振動数 In = 1; %正弦波の振幅

omega = 1.5; %正弦波の角周波数 dt = 0.01; %時間刻み

T = 60; %最終時間

xinit =[0;0]; %初期値セット

para = [xi,wn,In,omega]; %パラメータセット

%@SOE_FRで関数のポインタの受け渡し

[t1,xy] = RK2(@SOE_FR,[0,T],xinit,dt,para); %2次 Runge-Kutta

x = xy(1,:); xmin = min(x)-0.1; xmax = max(x)+0.1;

y = xy(2,:); ymin = min(y)-0.1; ymax = max(y)+0.1;

subplot(2,1,1); %2つのグラフを縦に並べて表示

plot(t1,x); %正弦波に対する応答の時間発展をプロット

axis([0 T xmin xmax]);

xlabel('t'); ylabel('x');

subplot(2,1,2); %2つのグラフを縦に並べて表示 plot(x,y); %状態空間プロット

axis([xmin xmax ymin ymax]);

xlabel('x'); ylabel('dx/dt');

周波数応答

正弦波入力

時間が十分経つと

€

x t ( ) ^{= G j} ( ) ^{ω sin} ( ^{ω t + tan}

⁻¹

^{G j} ( ) ^ω )

€

G j ( ) ω ^{= ω}

ⁿ

2

− ω

²

+ 2 j ζω

_n

ω + ω

_n²

€

d

²

x t ( )

dt

²

+ 2 ζω

_n²

^{dx t} ( )

dt + ω

_n²

^{x t} ( ) ^{= ω}

ⁿ²

^sin ( ) ^{ω t}

€

= 1

− ω

²

ω

_n²

+ 2 j ζ ω ω

_n

+ 1

ボード線図

€

ω

ω

_n

配布プログラム

•  伝達関数  TF_SOE.m

•  メイン部分 FOL_boad_main.m  

function [ y ] = TF_SOE(s,para)

y = para(2).^2./(s.^2 + 2.*para(1).*para(2).*s +…

para(2).^2.*ones(size(s)));

end

伝達関数 （ TF_SOE.m ）

パラメータ

出力 複素数s

Matlab は複素数が扱える

sがベクトルの場合に対応

€

G s ( ) ^{= ω}

ⁿ

2

s

²

+ 2 ζω

_n

^s + ω

_n²

メイン部分

SOE_boad_ma in.m

減衰係数を変えて、

周波数応答の変化を 観察

function main

xi = 0.05;%減衰係数 wn = 1;%固有振動数 j=sqrt(-1);%虚数の設定

para = [xi,wn]; %パラメータのセット

omega = 10.^(-1:0.05:1); %角周波数のセット 

G = TF_SOE(j.*omega,para); %周波数伝達関数の計算 gain = 20.*log10(abs(G)); %ゲインの計算

kaku = angle(G); %位相の計算

subplot(2,1,1); %2つのグラフを縦に並べて表示

semilogx(omega,gain); %x軸を対数にしてグラフ表示 xlabel('omega'); ylabel('dB'); %軸のラベル

subplot(2,1,2); %2つのグラフを縦に並べて表示

semilogx(omega,kaku); %x軸を対数にしてグラフ表示 xlabel('omega'); ylabel('rad'); %軸のラベル

end

レポート

•  関数 RK2.m を改造し、 4 次の Runge-Kutta 法 による数値積分を実行する関数を作成せよ。

•  作成した 4 次の Runge-Kutta 法の関数を用

いて、 2 次振動系のステップ応答を数値計算

し、グラフを描け。

ドキュメント内 dx dt = f x,t ( ) t (ページ 32-47)