Processing math: 0%

2017年7月10日月曜日

pythonでネイピア数を計算する

今日はpythonでネイピア数を計算して見ます。

ネイピア数は指数関数a^xを微分しても、元のa^xになるようにaを選んだもので、その時のaeと記述しネイピア数と呼んでいます。式で表すと以下のようになります。

\begin{equation} \frac{d e^x}{d x} = e^x \end{equation}

 ネイピア数は上の微分の式から、微分の定義に従って求めてみると、以下の式で表すことができることが、わかります。

\begin{equation} e= \lim_{n \to \infty} \left( 1+\frac{1}{n}\right) ^n \end{equation}

 数値計算でnをだんだん大きくして求めていきその様子をグラフ化して見ます。

#ここからソース
import matplotlib.pyplot as plt
import numpy as np
x=np.linspace(0,100,1000)
y=(1+1/x)**x
plt.plot(x,y)
plt.grid()
plt.show()
print(y[-1])
view raw napier.py hosted with ❤ by GitHub
#ここまでソース


2017年7月4日火曜日

pythonで1次遅れのボード線図を描く その2


一般的な場合の1次遅れの伝達関数


前回、具体的な伝達関数を使って、1次遅れの伝達関数のボード線図を書いて見ました。
では、一般的な場合はどうかと言うことを考えて見たいと思います。

さて、一般的な場合というのはどういうことかというと、1次遅れの伝達関数をゲインKと時定数\tauというものを使って記述するのが、普通のようです。
以下のようになります。

\begin{equation} G(s)=\frac{K}{\tau s + 1} \end{equation}

周波数伝達関数


さて、この式からボード線図を描くのに必要なゲインと位相を求めましょう。
まずsj \omegaを代入します。そして分母を有理化します。

\begin{equation} G(j \omega)=\frac{K}{1 + \tau \omega j} \frac{1- \tau \omega j}{1-\tau \omega j} = \frac{K(1 - \tau \omega j)}{1+(\tau \omega)^2} \end{equation}

ゲインの計算


ゲインを求めると

\begin{equation} |G(j \omega)|=\sqrt{\frac{K^2 + (-K \tau \omega )^2}{\{1+(\tau \omega)^2\}^2}} =\sqrt{ \frac{K^2 \{1 + (\tau \omega)^2 \}}{\{1+(\tau \omega)^2\}^2} } =\frac{K \sqrt{1+ (\tau \omega)^2}}{ 1 + (\tau \omega)^2} \end{equation}

\begin{equation} \begin{split} 20 \log |G(j \omega)|&=20 \log \frac{K \sqrt{1+ (\tau \omega)^2}}{ 1 + (\tau \omega)^2}\\ &=20 \log K + 10 \log \{1 +  (\tau \omega)^2 \} - 20 \log \{1 +  (\tau \omega)^2 \}\\ &=20 \log K - 10 \log \{1 +  (\tau \omega)^2 \} \end{split} \end{equation}

簡単になりましたね。

位相の計算


位相を求めましょう。

\begin{equation} \angle |G(j \omega)| = \tan^{-1} ( -\tau \omega)= - \tan^{-1}  \tau \omega \end{equation}

数値計算


#ここからソース
import numpy as np
import matplotlib.pyplot as plt
tau = 2.0
K=5
omega=np.logspace(-2,2,100)
gain=20*np.log10(K) -10*np.log10(1+(tau*omega)**2)
phase=-np.arctan(tau*omega)*180/np.pi
plt.subplot(2,1,1)
plt.plot(omega, gain,'-')
plt.xscale('log')
plt.grid(which='major')
plt.grid(which='minor')
plt.ylabel('Mag[dB]')
plt.xlabel('Freq[rad/s]')
plt.subplot(2,1,2)
plt.plot(omega, phase, '-')
plt.xscale('log')
plt.grid(which='major')
plt.grid(which='minor')
plt.ylabel('phase[deg]')
plt.xlabel('Freq[rad/s]')
plt.show()
#ここまでソース