コンテンツへスキップ

研究などで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)「完成」のコードを変更しました

参考

Scilab/Xcosが必要になったので入れて使ってみます。動作の確認として、MATLABでやったのと同じ解析をしてみました。

Scilab/Xcosが必要になったので入れて使ってみます。windows10,Scilab6.0.2です。動作の確認として、MATLABでやったのと同じ解析をしてみました。

目次

  • 導入
    • scilabとは
    • Xcosとは
  • インストール
  • 使ってみる
    • ブロックを並べる
    • 設定
    • シミュレーション
    • 構文
      • polyとは
      • syslinとは
      • bodeとは
      • dampとは
  • MATLAB/SimlinkとScilab/Xocsの比較
    • MATLAB/Simlink
    • Scilab/Xcos

導入

Scilabとは

Scilabはフリーのオープンソフトウェアで、端的に言ってしまえばMATLABのフリー版です。様々な数値解析ができます。

詳しくは公式にて

MATLABに関してはこちらの記事をご覧ください。=>MATLAB使い方個人的まとめ

Xcos

動的な連続/離散時間システムのモデリングやシミュレーションができます。 MATLABにとってのSimlinkです。グラフィカルにブロックを並べることでモデリングでき、コンパイルやシミュレーションもまとめてできる優れものです。

X2C(https://www.lcm.at/virtuelle-entwicklung/x2c/)

インストール

公式ページ: Download Scilab 6.0.2 -scilab

言語設定でjapaneseを選択し、デフォルトのままfullバージョンでインストールします。

インストールページのスクショ

使ってみる

MATLABと同様に、コマンドウィンドウ、Xcosウィンドウ、SimNoteのウィンドウをそれぞれ使って解析を行います。アプリケーション/Xcosで起動するとXcosウィンドウとライブラリが現れます。

Scilabの画面
Xcos

ブロックを並べてブロック線図を作り、sin波を表示させてみます。

ブロックを並べる

まずは信号源、出力/表示などからsin波発生器、CSCOPE、CLOCKを見つけて右クリックかドラッグでXcos上に配置し、矢印をドラッグして配線します。 スコープは赤い矢印があるものを選びましょう。これはイベントシミュレーションを駆動するイベント信号になります。

ブロック図の作成

設定

それぞれのブロックで右クリックすると項目が現れて、その中からブロックパラメータを開きます。 CFSCOPEを設定します。[参考1]

  • Color or mark vector: スコープで表示するグラフ線の色またはマークの連続によるグラフを指定する。入力信号数にかかわらず8セットの数字が必要。数字が正のとき、連続線で色が指定でき、負のとき、マーク形状が指定できる。(plot2dコマンドに準ずる)
  • Output window number: グラフ・ウィンドウ番号の指定。
  • Output window position: グラフ・ウィンドウの表示位置指定。左上の位置をx,yのドットで示す。
  • Output window size: グラフ・ウィンドウの表示サイズをドット数で指定。
  • Ymin: グラフ・ウィンドウのY軸最小値
  • Ymax: グラフ・ウィンドウのY軸最大値
  • Refresh period: グラフ・ウィンドウのX軸(時間)最大値 シミュレーション設定時間がこの値より大きい場合、グラフをクリアし再度続きの時間から描画する。基本的使い方では、シミュレーション時間と一致させる。
  • Buffer size: 入力データを一時的に蓄えておくメモリ領域のサイズを指定。グラフ表示速度に大きく影響するので、変化を観察したい場合以外は、少なくともRefresh period以上に設定する。最小は2。
  • Accept herited events: (inherited?) 継承の無効(0)・有効(1)を指定。1にするとイベント・ポート赤が削除される。アクティベーションが上流で行われる場合、継承が可能。

SCOPEの積算時間を10に、最小値~最大値を-10~10にしました。 クロックとsin発生器のパラメータも設定します。初期化時間を0.0に、ゲイン(振幅)5にします。

シミュレーションの設定

シミュレーション

シミュレーションを開始します。シミュレーション/パラメータ設定から、[参考1]

  • Final integration time: シミュレーション終了時間[単位は任意]
  • Real time scaling: 正の実数nを入力すると、シミュレーション設定単位時間tがn×t秒としてシミュレーションが実行される。n=1とした場合、シミュレーションの設定単位時間が1秒として解析表示される。
  • Integrator absolute tolerance: ソルバーの絶対誤差公差
  • Integrator relative tolerance: ソルバーの相対誤差公差
  • Tolerance on time: ODE(常微分方程式)ソルバーの最小時間公差
  • Max integration time interval: ソルバー呼び出しの最大時間間隔。"too many calls"というエラーメッセージが出た場合、値を小さくする。(解析精度を落とす。)
  • Solver 0(CVODE) - 100(IDA): 数値ソルバーの選択。
  • maximum step size: ソルバーがとる最大時間ステップ

最終積算時間を1.0E1.0、すなわち10sにしました。 AEBならA*10^Bになります。

document/test_0620/image.scgはこのようになりました。指定した通りに動いていることが確認できます。

document/ test_0620/image.scg の出力

構文

ビギナーのための科学技術計算ソフト(Scilab)の使い方講座 3. 付録: Scilab言語の基本構文 -SimCircuit

Scilab コンソールの使い方の基本

「;」を後ろにつけると結果を表示しないみたいですね。

polyとは

[参考2]によれば、ploy関数は以下のように書きます。

p=poly(実数もしくは行列,"変数")
-->a=[-3 -4];
-->p=poly(a,"x")
 p  = 

            2
   12 +7x +x 

コンソールに上のように入力すると、12+7x+x^2を返します。((x+3)(x+4)=0, のときx=-3,-4) pを多項式変数として扱う事ができます。

-->s=poly(0,'s')
s =

    s
polyをコンソールに入力したときの出力

syslinとは

[参考3]によれば、syslin関数は以下のように書きます。

q=syslin(連続時間システムc or離散時間システムd,多項式)
--> s=poly(0,'s');

--> P=syslin('c',1/(s+1))
 P  = 

           
     1     
   ------  
           
   1 + s   

bodeとは

[参考4]によれば、syslin関数は以下のように書きます。

bode(syslinで定義された線形システムの多項式,x軸の最小値,x軸の最大値,x軸の刻み,['コメント'])

コメントにタイトルなどをかきます。 こだわらないなら関数のみ与えればOKです。

dampとは

[参考5]によれば、syslin関数は以下のように書きます。

[wn,z] = damp(連続時間のsys関数)

MATLAB/SimlinkとScilab/Xcosの比較

某授業のレポート1でやった問題の伝達関数を使って確認していきます。

問題: オープンループの伝達関数L(s)が以下の時の、回ループと閉ループでのボーデ線図をそれぞれ表示します。

$$L(s)=\frac{0.1}{s(s+0.5)(s+1)}$$

MATLAB

  • オープンループ

$$L(s)=\frac{0.1}{s^3+1.5s^2+0.5s}$$

MATLABオープンループでのボーデ線図

位相余裕:59° 利得余裕:17.5dB

  • フィードバック $$L(s)=\frac{0.1}{s^3 + 1.5s^2 + 0.5s + 0.1}$$

A. 周波数応答

MATLABフィードバックでのボーデ線図

B. ステップ応答

設定がうまくいかず、ぎこちないですがこんな感じ。最終値は1です。

MATLABフィードバックでのステップ応答

Scilab

  • オープンループ

複数行のコマンドを何度も使用する場合、アプリケーション/SciNoteに書くとスクリプトのようにまとめて実行できます。 見やすさのために出力なしで実行しました。

SciNoteを使用している様子

A. 周波数応答

まずはコマンドで実行してみます。

コード

s=poly(0,'s');
p=syslin('c',0.1/(s*(s+0.5)*(s+1)))
bode(p)
show_margins(p)
[gm,fp]=g_margin(p) disp("ゲイン余裕[dB]") disp(gm) disp("位相交差角周波数[Hz]") disp(fp) 
[pm,fg]=p_margin(p) disp("位相余裕[°]") disp(pm) disp("ゲイン交差角周波数[Hz]") disp(fg) 

出力

ゲイン余裕[dB]

   17.501225

 位相交差角周波数[Hz]

   0.1125395

 位相余裕[°]

   59.289834

 ゲイン交点角周波数[Hz]

   0.0293667

ゲイン余裕:位相が-180°と交差する周波数(位相交差角周波数)においてゲインが0dBより何dB小さいか

位相余裕:ゲインが0dBと交差する周波数(ゲイン交差角周波数)

オープンループボーデ線図は

Scilab オープンループでのボーデ線図

MATLABと一致してますね。

B. ステップ応答

Scilab オープンループでのステップ応答
  • フィードバック

A. 周波数応答

コマンドからやります。フィードバックした場合の伝達関数G(s)は以下のようになるので

$$G(s)=\frac{L(s)}{1+L(s)}$$

コード

s=poly(0,'s');
p=syslin('c',0.1/(s*(s+0.5)*(s+1)))
g=p/(1+p)
disp(g)
bode(g)
show_margins(g)
[gm,fp]=g_margin(g) disp("ゲイン余裕[dB]") disp(gm) disp("位相交差角周波数[Hz]") disp(fp) 
[pm,fg]=p_margin(g) disp("位相余裕[°]") disp(pm) disp("ゲイン交差角周波数[Hz]") disp(fg) 
[wn,z] = damp(g); disp("固有周波数[Hz]") disp(wn) disp("減衰係数") disp(z) 

出力

                           
            0.1            
   ----------------------  
                    2   3  
   0.1 + 0.5s + 1.5s + s   

 ゲイン余裕[dB]

   16.258267

 位相交差角周波数[Hz]

   0.1125395

 位相余裕[°]

   114.71226

 ゲイン交差角周波数[Hz]

   0.0313482

 固有周波数[Hz]

   0.2964606
   0.2964606
   1.1378001

 減衰係数

   0.6108736
   0.6108736
   1.

同じ曲線を描いていることを確認しました。

Scilab フィードバックでのボーデ線図

B. ステップ応答

最終値1で、100秒まで計算しています。

Scilab フィードバックでのステップ応答

コードは[参考6]を参考にさせて頂きました。

雑記

  • Xcosのライブラリで固まったらEscキーを押しましょう!ドラッグより右クリックの方がいいです。
  • 周波数応答を調べるとき、関数に伝達関数を与えています。なのでXcosでブロック線図を描こうと思ったらブロックを並べて伝達関数を計算させて、それを取得して関数に入れる、という手順を踏むと思うので回りくどいです。ここでやった例はすぐ伝達関数を計算できるので、コマンドでやった方が早いと思いました。
  • 何気に変数ブラウザ便利ですね。

参考

1.ビギナーのための科学技術計算ソフト(Scilab)の使い方講座 2.3 Xcosによるシミュレーション・ツールとしての使い方 (GUI型ブロック線図入力モジュール) -SimCircuit

2. Scilab関数の説明:polyとは

3. Scilab関数の説明:syslinとは

4. Scilabでのボード線図の描き方

5. damp -scilab

6. scilabで制御系設計