図12.11のプロット

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

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

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

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

#図12.11のプロット

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

#図12.11(a)

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

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

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

#図12.11(b)

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

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

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

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

#図12.11(c)

#応答の初期値を与える
u3 = A*np.sin(omega3*t) #ω=10の正弦波

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

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

#図12.11(c)のプロット
plt.subplot(2,2,3) #複数の図を並べるためのコマンド．2行2列の2行1列目という意味
plt.plot(t3, y3, label="output") #周波数応答をプロット
plt.xlim([40,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("output") #縦軸のラベル表示
#plt.legend() #凡例の表示

#図12.11(d)

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

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

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