図12.12のプロット

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

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

#G(s)の分子・分母多項式を与える
num = [0, 0, 1] #分子多項式
den = [1, 0.3, 1] #分母多項式

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

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

#図12.12(a)

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

#正弦波（入力）の計算
u1 = A*np.sin(omega1*t) #ω=0.1の正弦波  
   
#周波数応答（出力）を求める
y1, t1, x1 = matlab.lsim(sys, u1, t, x0) #周波数応答

#図12.12(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.5,1.5]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("input and output") #縦軸のラベル表示
plt.legend() #凡例の表示

#図12.12(b)

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

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

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

#周波数応答（出力）を求める
y2, t2, x2 = matlab.lsim(sys, u2, t, x0) #周波数応答

#図12.12(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.5,1.5]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time   t[s]") #横軸のラベル表示
#plt.ylabel("input and output") #縦軸のラベル表示
plt.legend() #凡例の表示

#図12.12(c)

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

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

#正弦波（入力）の計算
u3 = A*np.sin(omega3*t) #ω=1の正弦波   

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

#図12.12(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([-4,4]) #縦軸の範囲の指定
plt.yticks([-4,-2,0,2,4]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("input and output") #縦軸のラベル表示
plt.legend() #凡例の表示

#図12.12(d)

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

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

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

#図12.12(d)のプロット
plt.subplot(2, 2, 4) #複数の図を並べるためのコマンド．2行2列の2行2列目という意味
plt.plot(t4, y4, label="output") #周波数応答をプロット
plt.xlim([0,50]) #横軸（時間軸）の範囲の指定
plt.ylim([-1.0,1.0]) #縦軸の範囲の指定
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() #凡例の表示
plt.show() #グラフの表示