コンテンツへスキップ

研究などでpythonを使って計算したりグラフ描画したりしたい人のための基礎事項まとめ本」でした。

「研究などでpythonを使って計算したりグラフ描画したりしたい人のための基礎事項まとめ本」でした。

科学技術計算のためのPython入門――開発基礎,必須ライブラリ,高速化 技術評論社 中久喜 健司 (著)

期間: 2019年9月4週目~10月1週目くらい

内容(個人的まとめ)

  1. pythonプログラミングの例
  2. IPython,Spyderなど
  3. pythonの基礎事項
  4. NumPy
  5. pandas
  6. SciPy
  7. Matplotlib
  8. 高速化

特徴

科学技術計算でPythonを使用したいと考えている方に向けた本になっていました。特定の何かをしようというものではなく、例が豊富で実行しながらpythonでできることを一通り網羅的に紹介して貰える本です。

感想

(参考として)読む前:pythonのフロー制御とモジュールを使って簡単な計算を自動化するプログラムが作成できる程度の知識がありました。一方で、Open CVなどその時々に必要なモジュールがないかググり、Qiita等に書いてある範囲内で使用することしかできませんでした。

全体

読む以前から多少の知識がある状態で、基本事項を改めて確認し直したい私にはとてもちょうどいい本でした。もしかしたらpythonを書いたことのない方には難しいかもしれませんが、例が豊富なのでひとまず読むことはできるかと思います。

公式リファレンスを読むのはまだ重い...という私と同じような方におすすめです。

1. pythonプログラミングの例

プログラムを書き、テストし、高速化して利用するという流れを一通り見せてくれようとしている章でした。作業フローを教えてくれる本はあまりであったことがなかったので貴重です。

例としてロケットシミュレーターを挙げており、公開されているコードを落としてこれるので実行しながら読み進められるはずです。が!地図をプロットするmatplotlib.basemapがエラー吐きまくりだったのでコードの実行はあきらめました。

2. IPython,Spyderなど

私はJupyter Notebookで普段書いています。Jupyter Notebookで書くとグラフや画像がその場で表示されるしTab補完もできるのでで重宝しています。

IPythonとはどんなものか、Spyderとは何かといった、知らなくてもどうにかなるが知っていると便利、でもわざわざ調べないことが書いてあり有難かったです。

3. pythonの基礎事項

一応知っていた内容が多い部分でした。「インデキシングとスライシングっていう名前の上での分類があったんだー」というような、教科書改めて読むと面白いという感じの内容です。

4. NumPy

NumPyいつもimportするけどnp.piしか使わないじゃん、と思っていた俺を殴りたい。多次元配列ndarrayとはどのようなものか、何がいいのかということを、例を実行しながらイメージできました。

5. pandas

よく聞くpandasとNumPyの違いとは何か、どんな時にpandasを使えばいいのか、なんとなく分かった気になれます。主にデータ型Series,DataFrame,Panelsに関して書かれています。

ただし、Panelsに関しては、「将来的になくなるよ、使わない方がいいよ」と警告してもらったので(下記)読み飛ばしました。本の悪い部分は情報の更新がわからないことなので仕方がないですね。

Panel is deprecated and will be removed in a future version.
The recommended way to represent these types of 3-dimensional data are with a MultiIndex on a DataFrame, via the Panel.to_frame() method
Alternatively, you can use the xarray package http://xarray.pydata.org/en/stable/.

6. SciPy

SciPyの項目では、以下の項目について例を挙げていました。

  • 統計関数
  • 離散フーリエ解析
  • ボーデ線図
  • データの内挿(グラフに近似線を追加するなど)
  • デジタル信号フィルタの設計
  • 行列の分解

7. Matplotlib

使いたいときにググればいいと思い、読み飛ばしました。

8. 高速化

自分が今必要としていない部分だと思いましたので、読み飛ばしました。

次に読みたい

関連記事

CSVファイルを指定して極座標系でグラフを描き、画像に出力するまでをやります。

ファイルを指定して極座標系でグラフを描き、画像に出力するまでをやります。

目次

  • 作業環境
  • CSV
  • 極座標系
  • ファイルのインポート

作業環境

  • anaconda 1.7.2
  • python 3.7.3
  • vscodeとexcel

jupyter notebookというIPythonを使用します。

jupyter notebook

jupyter notebook はサーバを立ててwebブラウザ上で文章を表示させるものです。 テキストエディタと同じ感覚で使用でき、IPythonの機能であるコードの解説の表示や図を含む実行結果をまとめて表示できます。

図を描くときは、実行してすぐ図を確認できるjupyterが楽です。anacondaをインストールするとまとめてインストールしてくれます。

CSV

csvって何

カンマで区切られた数値です。

VS code上で見やすく

csvはカンマで区切られているだけのデータなので見にくいです。そこで私が使っているvscode拡張機能を紹介します。

Rainbow CSV

Rainbow CSV

Excel VIewer

Excel viewer

ファイルを右クリックかcommand palletでpreviewが見れます。

極座標系の表示

=>matplotlib.pyplot.polar -matplotlib

import matplotlib.pyplot as plt
import numpy as np

r = 1
theta = 45
theta_rad = np.deg2rad(theta)# 角度をradに変換

ax=plt.subplot(projection="polar") #軸を表示
ax.scatter(theta_rad, r, c="red", s =15)

theta = np.arange(0.0, 4*2*np.pi, 0.01) #角度
r = 0.5*theta #半径
plt.polar(theta,r) # 極座標グラフのプロット

plt.savefig('PolarCoordinate.png')
plt.show()

subplotについて

fig.add_subplot(行,列,場所)を表します。(1,2,1)では列が(つまり縦方向に)半分になったグラフのうち1つめ(上側)が1つ表示されるらしいです。

import matplotlib.pyplot as plt

fig = plt.figure()
fig.add_subplot(1,2,1)
plt.show()
fig.add_subplot(1,2,1)
ax1 = fig.add_subplot(2,2,1)
ax2 = fig.add_subplot(2,2,2)
ax3 = fig.add_subplot(2,2,3)

ファイルのインポート

使っているデータは私が他の目的で作ったデータを適当に弄った意味のない数値列です。

import pandas as pd

df = pd.read_csv('fact_x2.csv')
print(df)
     turn[] length[mm]   total [mm]  radius[mm]    deg[deg]    \
0         0       0.00         0.00      0.0000    0.000000     NaN   
1         1      11.45        11.45      7.4000   88.653605     NaN   
2         1       0.45        11.90      7.4000   92.137808     NaN   
3         1      11.45        23.35      7.4000  180.791412     NaN   

~割愛~

117      17      16.10       854.45      9.4496   52.276096     NaN   
118      17       1.15       855.60      9.4496   59.248893     NaN   
119      17      16.10       871.70      9.4496  156.868057     NaN   

[120 rows x 10 columns]
import pandas as pd

df = pd.read_csv('fact_x2.csv')
print(df['turn[]'])
0       0
1       1
2       1
3       1

~割愛~

117    17
118    17
119    17
Name: turn[], Length: 120, dtype: int64

​指定できるのは行(Colum)。 横に伸びてる方です。

与える角度はdegじゃなくてradですので注意。

完成

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

df = pd.read_csv('fact_x2.csv')
r = df['Radius[mm]']
deg = df['Deg[deg]']
radian = np.deg2rad(deg)

fig = plt.figure(figsize=(10, 10))#新しい図面
ax = fig.add_subplot(111, projection='polar')#軸の表示
ax.plot(radian,r, '-o')
i=2
while i < 122:
    ax.plot([radian[i],radian[i-1]],[r[i], r[i-1]], 'or-')
    i +=2
print(np.pi)
plt.savefig('PolarCoordinate.png')
plt.show()

関数化

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

"""
CSVファイルから極座標系のグラフをプロットする
入力:(CSVファイル名, 出力画像ファイル名,半径、角度)
出力:出力画像ファイル名.png
"""
def polar(in_file,out_file,param_r,param_theta):
    df = pd.read_csv(in_file) #ファイル読み込み
    r = df[param_r]
    deg = df[param_theta]
    radian = np.deg2rad(deg) #deg をradに変換
    Len=len(r)#配列の長さを取得

    fig = plt.figure(figsize=(10, 10))#新しい図面
    ax = fig.add_subplot(111, projection='polar')#軸の表示
    ax.plot(radian,r, '-o') 
    i=2
    while i < Len: #1つ飛ばしどうしを赤線で結ぶ
        ax.plot([radian[i],radian[i-1]],[r[i], r[i-1]], 'xr-')
        i +=2
    plt.savefig(out_file)
    plt.show()

polar('test.csv','Polar_test.png','Radius[mm]','Deg[deg]')

雑記

進捗がない!!!!!!!!!!!
[追記](2019/9/19)「完成」のコードを変更しました

参考

pythonでグラフを表示するまでの流れを紹介します。グラフは自由振動、減衰振動、強制振動、倒立振り子(倒立振子、縦棒制御)の制御っぽいことをして時間変化をグラフで表示しました。おまけでCR回路の電圧の充電具合をプロットしました。 グラフを表示するのはCとかでも簡単にできますが、機械学習やってみたいのでpythonに触れる機会にします。

目次

  • 目的
  • グラフの様子
  • pythonとは
  • 環境とインストール
  • Jupyter Notebookの使い方
  • 振動の基礎とコード
  • 倒立振り子の運動方程式
  • 雑記
  • 参考

目的

思いのままのグラフを表示します。多項式だけではなく行列やベクトル、微分方程式も扱えるようになりたいと思います。

グラフの様子

今回は2つの似たような数式モデルをプロットしました。

図1:倒立振り子
  • 倒立振り子(縦棒制御)

手の上で箒のバランスを取るやつです。振動の中でも強制振動です。ニュートンの運動方程式を立てると、次のようになります。振り子は2l[m]、棒はm[kg]、台車はM[kg]です。このとき、角度Θを十分に小さいとして 近似してかなり簡単な式にまとめてあります。導出過程はページ最後に書きました。

$$\dot \theta_1(t)=\theta_2 \\ \dot \theta_2(t) = \frac{3}{l(4M+m)}(-L\theta_2+{g(M+m)-K}\theta_1)$$

簡単のためにl=1[m],M=2/3[kg],m=1/3[kg]を選ぶと、

$$\dot \theta_1(t)=\theta_2 \\ \dot \theta_2(t) = -L\theta_2+(g-K)\theta_1$$

この式で常微分方程式を解くと下記のような美しいグラフに!!コードが汚いのは御愛嬌ということで。python初めてだったのでほぼパクリです。減衰比については振動の基礎の項目に軽く書いています。

図2:倒立振り子の角度制御
from sympy import *      # sympyライブラリから全ての機能をimoprt
import numpy as np       #ベクトルや行列を表すのに使うライブラリnumpy
from scipy.integrate import odeint #scipyに含まれる常微分方程式を解く関数
import matplotlib.pyplot as plt #グラフを書くためのライブラリ
%matplotlib inline     #jupyter notebook内で表示させるため

"""
2階微分方程式 倒立振子 
"""
#シンボル定義
x=Symbol('x')                  # 文字'x'を変数xとして定義
t=Symbol('t')                  # 文字 't'を変数tとして定義

g = 9.8     #重力加速度
La = 0.89   #減衰比0.7のグラフのフィードバックゲイン
Ka = 10.2  
Lb = 1.0    #減衰比1.0のグラフのフィードバックゲイン
Kb = 10.05  
Lc = 0.0    #減衰比0.0のグラフのフィードバックゲイン
Kc = 10.0   
Ld = 1.0    #減衰比0.3のグラフのフィードバックゲイン
Kd = 12.57  

def func1(x, t, K):
    x1,x2=x
    dxdt=[x2,(-La*x[1])+((g-Ka)*x[0])]
    return dxdt
#同じことを複数書いているので省略

x0 = [10.0,2.0] # 初期条件 [x(0), dx(0)/dt]を表す

t=np.linspace(0,30,1000) # 時間 0から10までを101等分する
sol1=odeint(func1, x0, t, args=(k,)) # 微分方程式を解く 
#同じことを複数書いているので省略

#グラフ表示
plt.plot(t, sol3[:,0], linewidth=1,label='ζ=0.0') #  角度を図示。
plt.xlabel('time t[s]', fontsize=18)
plt.ylabel('angle θ[rad]', fontsize=18)
plt.legend(loc='upper right')
plt.savefig('huriko1.jpg')
plt.show()
  • LCR回路

そしてもう一つの式がこれです。回路のキャパシタはダンパー、インダクタはばねと同じです。

$$ E = v_R + v_C + v_L \\ E = L\frac{di(t)}{dt}+Ri(t)+\frac{1}{C}\int i(t)dt $$

これをキャパシタ電圧でまとめると

$$ i = C\frac{d}{dt}v_C \\ \frac{d^2}{dt^t} v_c = -\frac{R}{C} \frac{d}{dt}v_c -\frac{1}{LC}v_c + \frac{E}{LC}$$

これは強制振動の式と同じ形なので、同じようなグラフがかけます。

図3:LRC回路のVcの時間応答

python とは

今回のpythonは読みやすく、それでいて効率もよいコードをなるべく簡単に書けるように作られたプログラミング言語です。言語自体も書きやすいように作られていますが、ライブラリやワーススペースが豊富で、だいたい何でも出来るみたいです。最近ではAIや機械学習が人気ですが、そういう事をするためのライブラリが豊そろっているみたいで、最近活躍している言語です。 2系と3系は全然違うので注意が必要です。

ライブラリというのは、

ライブラリとは、複数のプログラムが共通して利用するコードをまとめたファイルである。 ライブラリ自体は単独で実行することはできない。 from I-8-8. ライブラリとプログラムの関係 | 日本OSS推進フォーラム

というものです。必要なときに呼び出すと自分の仕事だけやって帰って行きます。

ちなみに、初学者へ。世の中のホームページやアプリの画面表示は別な言語で書いています。pythonとかCとかは計算処理をやってくれるものなので、これだけやっても画面表示は出来ません。注意してください。勉強しはじめたころの私はこのあたりよく分かって無かったために時間が掛かってしまったので一応書いておきました。

環境とインストール

環境はAnacondaでまとめてインストールしてJupyter Notebookで書きます。 python公式でインストールしても良いんですが、他にもいろいろまとめてインストール出来るらしいDownloads - Anacondaでインストールした方が楽みたいです。デフォルトのままで良いので指示に従うだけです。

Jupyter Notebookの使い方

インストールしたらJupyterNotebookを起動します。するとコマンドコマンドプロンプト(っぽいなにか)が出てきますので、そこに書かれているURLをブラウザで開いてください。

ユーザーのフォルダの中身が表示されますので、好きなのを選んで右上のnewからpython3をクリック。すると、プログラムを書くためのセルが出てきます。

今数字が書いてある色塗りの部分を押すとコマンドモードになり数式やらがかけます。In[]やOut[]と書いている部分はエディットモードで、コピーやらが出来ます。エディットモードでHを押すとヘルプが出ます。

コマンドモードで上のcodeの部分をmorkdownにすると、morkdownで文章を書けます。 morkdownは私がこのブログの下書きやメモに普段使いしているマークアップ言語(HTMLの仲間)で、# ABC と書くとタイトル的な大きな文字、## ABCと書くと1段階小さい文字、なにもつけないと普通の文字、- と書くと箇条書き、というようにかなり使いやすい言語です。

pythonの記法はすっ飛ばします。numpyというライブラリを読みこんでベクトル(1次元配列)や行列(2次元配列)を表せます。 スライシング。bool。

コマンドプロンプトでサーバーをlocalhostで立てて、ブラウザで表示しているんじゃないかと思います(このあたりよく解ってませんが。たぶん。。なので、ブラウザで開いている間にJupyterNotebookを消すとエラーをはきます。

発生したエラーについて

  • 保存した画像が真っ白になっている。jupyterにはグラフが出力されている

Python matplotlibのグラフが保存できない(windows+anaconda) teratail

振動の基礎とコード

図4:線形振動系の力学モデル

機械屋さんにとってのかなり基本的なモデル(らしい)で考えます。ニュートンの運動方程式を立てます。

$$ f(t)= m\frac{d^2x(t)}{dt^2}+c\frac{dx(t)}{dt}+kx(t) $$

これは強制振動を表しています。

振動は、自由振動、減衰振動、強制振動の3つに分けられます。自由振動は正弦波などのように振動し続けます。減衰振動は自由振動が次第に減衰していきます。強制振動は減衰振動に外力を加えている場合です。 日常生活では大概空気抵抗などで減衰するので減衰振動がイメージしやすいかと思います。一方、制御するときはモータが動いているので強制振動です。

外力をf(t)と書いていますから、f(t)=0なら減衰振動、f(t)の項が有れば強制振動になります。

上の式を作為的に以下のように変形します。

$$ f(t)= \frac{d^2x(t)}{dt^2}+2\zeta\omega_n\frac{dx(t)}{dt}+\omega_n^2 \ \\減衰比\zeta = \frac{c}{2\sqrt{km}}, 固有振動数\omega_n = \sqrt{\frac{k}{m}} $$

減衰比ζは減衰具合、固有振動数ω_nは振動時の周波数に対応しています。これが減衰振動の基本的なモデルです。ζ=0のとき全く減衰しないのでsin波のような自由振動をし、ζ=1の時過不足無く(クッションに沈むように)目標値に到達します。実際には誤差や外乱が入るので制御系を作るときには0.7くらいにするそうです。

今回はpythonで解きたいので、ここから変形していきます。本当はいろんな方法が有るみたいなんですが、そのうちやると言うことにして、今回は参考本と同じ方法で解いていきたいと思います。

常微分方程式を解きます。二次方程式は大変なので、以下のように書きかえます。

$$\frac{dx_1(t)}{dt}=x_2 \\ \frac{dx_2(t)}{dt} = -\frac{c}{m}x_2-\frac{k}{m}x_1 +f(t)$$

これを行列を用いたベクトル表現にすると、

$$ \frac{d}{dt} \left[ \begin{array}{rrr} x_1 \\ x_2 \end{array} \right] = \left[ \begin{array}{rrr} 0 \ & 1 \ \\ -\frac{k}{m} & -\frac{c}{m} \\ \end{array} \right] \left[ \begin{array}{rrr} x_1 \\ x_2 \end{array} \right] + \left[ \begin{array}{rrr} 0 \\ 1 \end{array} \right] f(t)$$

この式のように変数の位置x(t)を配列としてpythonで定義すると、以下のようになります。

m = 1.0 #質量
c = 0.5 #粘性減衰
k = 1.0 #ばね定数

def func(x, t, k):
    x1,x2=x
    dxdt=[x2,-c/m*x[1]-k/m*x[0]+sin(t/2)]
    return dxdt

x0 = [0.0,0.0] # 初期条件 [x(0), dx(0)/dt]を表す.今は初速0.

t=np.linspace(0,50,100000) # 時間 0から10までを101等分する
sol=odeint(func, x0, t, args=(k,)) # 数値的に微分方程式を解き, x(t)とx'(t)をsolのリストの[:,0]、[:,1]成分へ格納する。

以上のような感じで、パラメータをいじることで自由振動、減衰振動、強制振動がそれぞれ再現出来ました。

  • 自由振動(図5) 初期x=0.0, 初期x' = 1.0, m = 1.0, c = 0.0, k = 1.0
  • 減衰振動(図6) 初期x=0.0, 初期x' = 1.0, m = 1.0, c = 0.5, k = 1.0
  • 強制振動(図7) 初期x=0.0, 初期x' = 0.0, m = 0.0, c = 0.5, k = 1.0 f(t) = sin(t/2)
  • 強制振動(図8) 初期x=0.0, 初期x' = 0.0, m = 1.0, c = 0.5, k = 1.0 f(t) = 1
図5:自由振動
図6:減衰振動
図7:強制振動
図8:強制振動

倒立振り子の運動方程式

再度掲載

台車と振り子を図に示します。 粘性摩擦係数をB、 軸の粘性摩擦係数をCとおきました。ここから、下の4つの式が立てられます。

振り子について

$$J\ddot \theta = Vlsin\theta - Hlcos \theta - C\dot \theta …(1)\\ (J=\frac{ml^2}{3})\\ m\frac{d^2}{dt^2}(x + lsin\theta)= H…(2) \\ m\frac{d^2}{dt^2}(lcos\theta)=V - mg…(3)$$

台車について

$$M \ddot x = f(t) - B \dot x - H…(4)$$

簡単のためにC = 0, B = 0とします。ここからH,Vを消去します。

$$H = m \ddot x + ml \ddot \theta sin\theta - ml (\dot \theta)^2 sin\theta…(5) \\ V = -ml \ddot \theta sin \theta - ml(\dot \theta) cos \theta +mg…(6)$$

これを(1)(4)に代入してまとめると、

$$ J\ddot \theta= -ml^2 \ddot \theta - lmcos\theta \ddot x + mglsin\theta…(7) \\ (M+m) \ddot x = -ml\ddot \theta cos\theta + ml(\dot \theta)^2 sin\theta+f(t)…(8)$$

この式からxの2階微分を消去して、Θが十分小さいとして線形化すると

$$cos\theta\approx 1,sin\theta\approx\theta, (\theta)^2 \approx 0 \\ \ddot \theta l (4M+m)-3(M+m)g \theta = -3f(t)…(9)$$

f(t)に角度と各速度についてフィードバックさせます。

$$ f(t)= L \dot \theta + K \theta\ $$

簡単のために l = 1[m], M = 2/3[kg], m = 1/3[kg]とおくと、以下のようにかけます。プログラムで解けるように一階化します。

$$ \dot \theta_1 = \theta_2 \\ \dot \theta_2 = -L \theta_2 + (g-K)\theta_1$$

一応全体の全体をラプラス変換して伝達関数の式を求めます。(せっかく導出したから書きたい。)

$$L{\theta(t)}=\Theta(s),L{f(t)}=F(s)として\\ \frac{\Theta(s)}{F(s)}=\frac{-3}{(4M+m)ls^2+3ls-3(M+m)g+K}$$

雑記

プログラミング系のサイトは最近公開されたものを見るように気をつけるようになりました。機械系や電気系の、特に理論はかなり成熟しているのでふるいサイトでも何でも読みやすくわかりやすければ良いんですが、情報系はツールがどんどん更新されていくので最新のを見ないと妙なところで躓いてしまいます。 おかげで私のPCちゃんはパスとかアプリとか訳が分からない感じになっています。情報系の特にweb系のブロガーさんは情報収集がすごいので、最新のでも誰かしら先にやっていますので安心ですね。先人は偉大です。 それと、解析力学が難しいですね。もう少し時間をかけないとわからないです。ラグランジュ難しい。制御はやって行きたいので、のんびり本を読んでいきます。
追記(2019/3/20)参考リンクをリストにしました。Amazonリンクを増やしました。

参考