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

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

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

$$\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

#ここまでソース

pythonで１次遅れのボード線図を描く その２

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

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

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

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

周波数伝達関数

さて、この式からボード線図を描くのに必要なゲインと位相を求めましょう。
まず$s$に$j \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()

view raw bode_first_order.py hosted with ❤ by GitHub

#ここまでソース

anaconda3でインストールしたmatplotlibのエラーの回避

Ubuntu１６．０４LTS環境で
anaconda3を用いてpython環境を構築しました。

matplotlibを使ってグラフを描こうとしたらエラーが出て止まります。
以前は出なかったので、あれこれ悩んで解決したので、備忘録として書いておきます。

２０１７年６月末時点でanaconda3でpython環境を構築するとPyQt5がインストールされるみたいですが、そのことがエラーの原因になってると思われます。

インタラクティブシェルでmatplotlib.pyplotをインポートしようとすると、以下のエラーが出ます。

ModuleNotFoundError: No module named 'PyQt4'
 

解決法ですが

matplotlibの設定ファイルmatplotlibrcを書き換えます。デフォルトではPyQt4用のなっていると推察されましたので、その部分を書き換えます。

設定ファイルがどこにあるのかわからないのですが、インタラクティブシェルで以下のように入力すると場所がわかります。

>>> import matplotlib

>>> matplotlib.matplotlib_fname()

場所が判ったらmatplotlibrcを適当なテキストエディタで開いて

元は以下のようになっていると思うので、赤いところを

#### CONFIGURATION BEGINS HERE
# The default backend; one of GTK GTKAgg GTKCairo GTK3Agg GTK3Cairo
# CocoaAgg MacOSX Qt4Agg Qt5Agg TkAgg WX WXAgg Agg Cairo GDK PS PDF SVG
# Template.
# You can also deploy your own backend outside of matplotlib by
# referring to the module name (which must be in the PYTHONPATH) as
# 'module://my_backend'.
backend      : Qt4Agg
# If you are using the Qt4Agg backend, you can choose here
# to use the PyQt4 bindings or the newer PySide bindings to
# the underlying Qt4 toolkit.
backend.qt4 : PyQt4        # PyQt4 | PySide

以下の青いように書き換えます

#### CONFIGURATION BEGINS HERE
# The default backend; one of GTK GTKAgg GTKCairo GTK3Agg GTK3Cairo
# CocoaAgg MacOSX Qt4Agg Qt5Agg TkAgg WX WXAgg Agg Cairo GDK PS PDF SVG
# Template.
# You can also deploy your own backend outside of matplotlib by
# referring to the module name (which must be in the PYTHONPATH) as
# 'module://my_backend'.
backend      : Qt5Agg
# If you are using the Qt4Agg backend, you can choose here
# to use the PyQt4 bindings or the newer PySide bindings to
# the underlying Qt4 toolkit.
#backend.qt4 : PyQt4        # PyQt4 | PySide

以上で、正常に動作するようになりました。

LaTeX備忘録

$\LaTeX$ を使っていてちょくちょく検索することをまとめておこうと思います。

括弧

\left( \frac{1}{s+1} \right)

$$\left( \frac{1}{s+1} \right)$$

\left\{ \frac{1}{s+1} \right\}

$$\left\{ \frac{1}{s+1} \right\}$$

pythonで一次遅れのボード線図を描く

pythonを使って制御工学の問題を解いていこうと思っています。少しづつ更新していきますので、今続きがなくても、たまに確認してみてくださいね。

ボード線図を描く

今回は一次遅れ系のボード線図を描くことを考えてみます。

伝達関数

まず一次遅れの伝達関数の一例は次の式で表されます。

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

周波数伝達関数

 ここでこの伝達関数のボード線図を描く前にすることは$s=j\omega$ と置き換えをします。

$$\begin{equation} G(j\omega)=\frac{1}{1+j\omega} \end{equation}$$


この置き換えをした後の伝達関数を周波数伝達関数と呼びます。

すると伝達関数は$\omega$を変数とする、複素数と見ることができます。

ところでボード線図は、ゲイン曲線と位相曲線の二つに分けて描きます。

ゲイン線図は周波数伝達関数からゲイン$20\log |G(j\omega)|$を計算したものを横軸が対数目盛の片対数グラフとしたものです。
横軸が角周波数$\omega$で、縦軸はゲインです。

複素数の実部と虚部そして有理化

複素数からゲインを求めるためには複素数の実部と虚部を明らかにしておく必要があります。以下ではその過程を述べてみたいと思います。
一次遅れの伝達関数は分母が複素数の形をした分数の形をした関数でこのような姿の関数を有理関数と呼びます。
複素数の分母を複素数ではないものにすることを有理化と呼びます。以下で有理化の作業を実際にやってみます。

有理化には分母の複素数の共役複素数を分母と分子にかけます。同じものを分母と分子にかけますから。かける前とかけた後も関数自身は同じものです。

$$\begin{equation} \begin{split} G(j\omega)&=\frac{1}{1+j\omega}\frac{1-j\omega}{1-j\omega}\\ &=\frac{1-j\omega}{1+\omega^2} \end{split} \end{equation}$$

以上が有理化と言われる作業です。複素数であった分母からは複素数がなくなりました。
さて、話を戻して実部と虚部を求めてみましょう。

複素数を眺めてみたときに虚数単位のjが付いていない項が実部で、jが付いてる項のjの係数を虚部と言います。
したがって、今の場合は、$\frac{1}{1+\omega^2}$が実部で、$-\frac{\omega}{1+\omega^2}$が虚部と言うことができます。

 ここで今後の記述が大変なので実部と虚部を以下のように置き換えておきます。

$$\begin{equation} \begin{split} x=&\frac{1}{1+\omega^2}\\ y=-&\frac{\omega}{1+\omega^2} \end{split} \end{equation}$$

これで、ボード線図を描くための準備がほぼ整いました。
ゲインと位相を計算します。

 ゲインの計算

ゲインを計算するためには複素数の大きさを求めなければなりません。大きさは以下の式から計算します。
$$\begin{equation} |G(j\omega)|=\sqrt{x^2+y^2}=(x^2 + y^2)^{\frac{1}{2}} \end{equation}$$

ゲインは以上の対数をとって２０倍したものと定義されています。

$$\begin{equation} 20 \log |G(j\omega)|=20 \log (x^2 + y^2)^{\frac{1}{2}}=10 \log (x^2 + y^2) \end{equation}$$

$x$と$y$に元の式を代入して整理してみましょう。

$$\begin{equation} 10 \log (x^2 + y^2)= 10 \log \left\{\left(\frac{1}{1+\omega^2}\right)^2 + \left(-\frac{\omega}{1+\omega^2}\right)^2\right\} \end{equation}$$

さて、さらに式の整理を進めましょう

$$\begin{equation} \begin{split} &10 \log \left\{\left(\frac{1}{1+\omega^2}\right)^2 + \left(-\frac{\omega}{1+\omega^2}\right)^2\right\}=10 \log \left\{ \frac{1+\omega^2}{(1+\omega^2)^2}\right\}\\ &=10 \log \frac{1}{1+\omega^2} \end{split} \end{equation}$$

だいぶ簡単になりましたね。ここで対数法則を使って、分数を引き算にします。

$$\begin{equation} 10 \log \frac{1}{1+\omega^2} = 10 ( \log 1 - \log (1+ \omega^2) ) \end{equation}$$

$\log 1$は０なので、さらに簡単になって以下のようになりました。

$$\begin{equation} 10 ( \log 1- \log (1+\omega^2) ) = -10 \log(1+\omega^2) \end{equation}$$

位相の計算

ゲインの計算がひとまず終わったので位相の計算をします。

位相は先ほど定義した実部$x$、虚部$y$を使って説明すると、以下のように表されます。

$$\begin{equation} \tan^{-1} \frac{y}{x} \end{equation}$$

ここでもに元の式を代入して計算を進めてみましょう。

$$\begin{equation} \begin{split} \tan^{-1} \frac{y}{x}&=\tan^{-1} \frac{-\frac{\omega}{1+\omega^2}}{\frac{1}{1+\omega^2}}\\ &=\tan^{-1} (-\omega)\\ &=-\tan^{-1} \omega \end{split} \end{equation}$$

 数値計算！

#ここからソース

	import numpy as np
	import matplotlib.pyplot as plt
	omega=np.logspace(-2,2,100)
	gain=-10*np.log10(1+omega**2)
	phase=-np.arctan(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()

view raw bode.py hosted with ❤ by GitHub

#ここまでソース

さて、上のソースが今回ボード線図を描くのに使用したソースコードです。
本記事はpythonでボード線図を描ききるのが主眼です。

import numpy as np
import matplotlib.pyplot as plt

でpythonに数値配列用のライブラリ超有名なnumpyを読み込みます。使うときはnpとして使います。

次に、グラフを描くためのライブラリであるmatplotlib.pyplotを読み込みます。
私は長年、仕事にgnuplotを使ってきましたが、だいぶ慣れてきたので、グラフ描画をpyplotでやってもいいかなと思い始めています。


描きあがったボード線図が以下です。

やったことを残す

最近、記録しておかないことはすべて無と同じという認識を持っている。

ものづくりをしてものができると、満足感が得られる。

しかし、ものは場所を取るし、時間が経つと劣化して、いつかは壊れる。だから、そのうち破棄することになり無くなってしまう。

無くなってしまった瞬間、そのものに関した記憶も無くなってしまうように感じる。

ものづくりだけではなく
仕事や勉強、その他諸々のものに関しても、記録がなければ無と同じかもと考えている。

三角比

　三角比は直角三角形の辺の比です。
　直角三角形には

　斜辺

　隣辺

　対辺

　と名付けられた３つの辺が有ります。

　斜辺以外の辺は、直角以外の２つの角のどちらを角度の基準にするかによって名前が隣辺又は対辺と変化します

上のように印の付いた角が何処かで、直角三角形の辺の名前を変えて考えると三角比を考えるときに便利です。