図12.5のプロット

#計算とグラフプロットに必要なモジュールの読み込み
import numpy as np
from control import matlab
from matplotlib import pyplot as plt
from scipy import arange 

#正弦波のパラメータを与える
A = 1.0 #正弦波の振幅
omega1 = 0.1 #ω=0.1の場合の正弦波の角周波数
omega2 = 0.5 #ω=0.5の場合の正弦波の角周波数
omega3 = 1.0 #ω=1の場合の正弦波の角周波数  
omega4 = 10 #ω=10の場合の正弦波の角周波数

#伝達関数のの分子・分母多項式を与える
num = [0, 0, 1] #G分子多項式
den = [2, 10.2, 1] #分母多項式

#伝達関数表現を与える
sys = matlab.tf(num, den) #伝達関数表現           

#時間変数の定義
t = arange(0, 50, 0.01) #0から50まで0.001刻み

#図12.5(a)

#正弦波（入力）の計算
u1 = A*np.sin(omega1*t) #ω=0.1の正弦波 

#応答の初期値を与える
x0 = [0, 0] #x(0)=[0, 0]と与える（2次遅れ系なのでベクトルとなる）        

#周波数応答（出力）を求める
y1, t1, x1 = matlab.lsim(sys, u1, t, x0) #周波数応答

#図12.5(a)のプロット
plt.subplot(2, 2, 1) #複数の図を並べるためのコマンド．2行2列の1行1列目という意味
plt.plot(t1, u1, label="input") #入力をプロット
plt.plot(t1, y1, label="output") #周波数応答をプロット
plt.xlim([0,50]) #横軸（時間軸)の範囲の指定
plt.ylim([-1,1]) #縦軸の範囲の指定
plt.yticks([-1.0,-0.5,0,0.5,1.0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("input and output") #縦軸のラベル表示
plt.legend() #凡例の表示

#図12.5(b)

#正弦波（入力）の計算
u2 = A*np.sin(omega2*t) #ω=0.5の正弦波

#応答の初期値を与える
x0 = [0, 0] #x(0)=[0, 0]と与える（2次遅れ系なのでベクトルとなる）           

#正弦波（入力）の計算
y2, t2, x2 = matlab.lsim(sys, u2, t, x0) #周波数応答

#図12.5(b)のプロット
plt.subplot(2, 2, 2) #複数の図を並べるためのコマンド．2行2列の1行2列目という意味
plt.plot(t2, u2, label="input") #入力をプロット
plt.plot(t2, y2, label="output") #周波数応答をプロット
plt.xlim([0,50]) #横軸（時間軸)の範囲の指定
plt.ylim([-1,1]) #縦軸の範囲の指定
plt.yticks([-1.0,-0.5,0,0.5,1.0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time   t[s]") #横軸のラベル表示
#plt.ylabel("input and output") #縦軸のラベル表示
plt.legend() #凡例の表示

#図12.5(c)

#正弦波（入力）の計算
u3 = A*np.sin(omega3*t) #ω=1.0の正弦波
    
#応答の初期値を与える
x0 = [0, 0] #x(0)=[0, 0]と与える（2次遅れ系なのでベクトルとなる）        

#周波数応答（出力）を求める
y3, t3, x3 = matlab.lsim(sys, u3, t, x0) #周波数応答

#図12.5(c)のプロット
plt.subplot(2, 2, 3) #複数の図を並べるためのコマンド．2行2列の2行1列目という意味
plt.plot(t3, u3, label="input") #入力をプロット
plt.plot(t3, y3, label="output") #周波数応答をプロット
plt.xlim([0,50]) #横軸（時間軸)の範囲の指定
plt.ylim([-1,1]) #縦軸の範囲の指定
plt.yticks([-1.0,-0.5,0,0.5,1.0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("input and output") #縦軸のラベル表示
plt.legend() #凡例の表示

#図12.5(d)

#正弦波（入力）の計算
u4 = A*np.sin(omega4*t) #ω=10の正弦波

#周波数応答（出力）を求める
y4, t4, x4 = matlab.lsim(sys, u4, t, x0) #周波数応答

#図12.5(d)のプロット
plt.subplot(2, 2, 4) #複数の図を並べるためのコマンド．2行2列の2行2列目という意味
plt.plot(t4, y4, label="output") #周波数応答をプロット
plt.xlim([40,50]) #横軸（時間軸)の範囲の指定
plt.ylim([-0.01,0.01]) #縦軸の範囲の指定
plt.yticks([-0.010,-0.005,0,0.005,0.010]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time   t[s]") #横軸のラベル表示
#plt.ylabel("output") #縦軸のラベル表示
plt.legend() #凡例の表示

plt.show() #グラフの表示