図11.2のプロット

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

#1次遅れ系のパラメータと正弦波A sin ωtのパラメータを与える
T = 1 #T=1
K = 1 #K=1

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

#正弦波A sin ωtのパラメータを与える
A = 1 #正弦波の振幅
omega1 = 0.01 #ω=0.01の場合の正弦波の角周波数

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

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

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

#応答の初期値を与える
x0 = [0] #x(0)=0

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

#図11.2(a)のプロット
plt.plot(t1, u1, label="input") #入力をプロット
plt.plot(t1, y1, label="output") #周波数応答をプロット
plt.xlim([0,1000]) #横軸（時間軸）の範囲の指定
plt.ylim([-1,1]) #縦軸（出力）の範囲の指定
plt.xticks([0,200,400,600,800,1000]) #横軸の目盛りの値の設定
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() #グラフの表示

#正弦波A sin ωtのパラメータを与える
A = 1 #正弦波の振幅
omega2 = 10 #ω=10の場合の正弦波の角周波数

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

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

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

#図11.2(b)のプロット
plt.plot(t2, u2, label="input") #入力をプロット
plt.plot(t2, y2, label="output") #周波数応答をプロット
plt.xlim([0,10]) #横軸（時間軸）の範囲の指定
plt.ylim([-1,1]) #縦軸（出力）の範囲の指定
plt.xticks([0,2,4,6,8,10]) #横軸の目盛りの値の設定
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() #グラフの表示

#図11.2(c)の出力（5〜10[s]の拡大図）

plt.plot(t2, u2, label="input") #入力をプロット
plt.plot(t2, y2, label="output") #周波数応答をプロット
plt.xlim([5,10]) #横軸（時間軸）の範囲の指定
plt.ylim([-1,1]) #縦軸（出力）の範囲の指定
plt.xticks([5,6,7,8,9,10]) #横軸の目盛りの値の設定
plt.yticks([-1.0,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1.0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("input and output") #縦軸のラベル表示
plt.legend() #凡例の表示
plt.show() #グラフの表示