個人的にPythonが有難いのは以下の点です。

①行列計算や多次元配列の取扱いが簡単

②複素数計算も可能（Bode線図やNyquist線図を作る際に便利）

③グラフ描画が容易で、綺麗な図が簡単に書ける

④外部ファイルとの入出力が行いやすい。(csv、txtは勿論、xmlやJSON等なんでも取り扱いやすい)

⑤変数宣言しないでガンガン書ける（この書き方はプログラミング専門の方から言わせれば問題外だと思いますが、、）

⑥COM接続、OPC接続、OSI PIとの接続等、様々な外部リンクが簡単に行える

⑦最適化計算や非線形方程式、常微分方程式までモジュールが充実している。

⑧Jupyternotebookを使うことで、技術情報とコードの出力結果をまとめた資料作成が可能

⑨controlライブラリも使用できる。

Pythonでの数値計算の恩恵を受けるためにいくつか便利なライブラリを以下に紹介します。

・SciPy、NumPy、matpoltlib、scikit-learn：最適化計算や逆行列演算、グラフ描画等に便利です。

・PySide：GUIの作成に使えます。

ちなみに私はインストールで何度もはじかれて苦戦しました。

インストールについてはゆくゆく整理していきたいと思います。

※2016.05.08追記：Scipy、Numpy、Matplotlibのインストールについては、こちらがよくまとまっています。

配列（array）、転置（transpose）、積(.dot)

※Numpyのインストールが必要です。

例えば以下の行列の表現は、x=np.array([[a,b,c],[d,e,f],[g,h,i]])となります。

$$x=\begin{bmatrix} a & b & c\\\ d & e & f \\\ g & h & i \end{bmatrix}$$

例えば以下の行列の表現は、a=np.array([[1],[2],[3]])となります。

$$a=\begin{bmatrix} 1 \\\ 2 \\\ 3 \end{bmatrix}$$

また上記の行列を転置した以下の行列はa=np.transpose(a)のように書けます。

$$a=\begin{bmatrix} 1 & 2 & 3 \end{bmatrix}$$

行列同士の積は、例えば行列aと行列bの積であれば、a.dot(b)と書くだけです。

例えば上記のシンプルな行列\(a=\begin{bmatrix} 1 & 2 & 3 \end{bmatrix}\)と、その転置行列の積をPythonで計算する場合、非常に短いコードで済みます。

コード例 　　 　　 　　 　　出力例

このようにNumpyを使うと、転置行列の作成、行列の積が簡単にできます。

csvファイルの読み込み

numpyのloadtxtを使うと簡単に読み込むことができます。

import numpy as np
a = np.loadtxt("sampledata.csv",delimiter=",")

numpyのloadtxtを使うと簡単に読み込むことができます。

csvファイルを読み込んで、データの統計処理を行うのも簡単です。

最大値

print(np.amax(a))

最小値

print(np.amin(a))

算術平均

print(np.mean(a))

中央値

print(np.median(a))

標準偏差

print(np.std(a))

分散

print(np.var(a))

Excelファイルの読み込みについて

numpyのloadtxtを使うと簡単に読み込むことができます。

import numpy as np
      a = np.loadtxt("sampledata.csv",delimiter=",")

scipy.optimeze

curve_fit()、leastsq()等、回帰計算に役立つ関数が入っています。

scipy.integrate

odeint()で簡単に微分方程式を記載することができます。

ode()を使うと、set_integrator()で数値解法を変更することが可能です。

使い勝手の良い関数をいくつかご紹介します。

read_csv(ファイルの置き場所,encoding)：csvファイルの読み込み

dropna()：欠損値の取り扱いに便利です。欠損値を含む(how=any)、あるいは全ての値(how=all)が欠損値である行、あるいは列(axis=1)を除外することができます。逆に欠損値を埋める場合はfillna()を使うことができます。

unique()：ユニークな値を持ったリストを作成することができます。

2dplot出力

詳細はゆくゆくまとめていきたいと思いますが、matplotlibのplotはsubplotをまとめて一つの形式で出力できるため、非常に便利です。

例えば、NRTLで水-EthanolのXY、PXY、活量係数をまとめてプロットしたいという場合、以下のようにまとめて表示させることが可能です。

以下のように各subplotの配置位置を2dgridで記載し、各suubplotの詳細を記載していきます。

3dplot出力

Qiitaのブログに詳細記載しましたのでご参照ください。

例として、以下のように理想気体と実在気体のPVT図を重ね書きする方法をまとめています。時間のあつ時にもう少し詳細に書きたいと思います。

三角座標のプロット

ブログに詳細記載しましたのでご参照ください。

例えば、以下のVapor pressureのデータ（データはperryのhandbook参照）をテキストファイルに保存しておき、各係数を必要に応じて読みだすことを考えます。

$$ln(P)=C1+C2/T+C3ln(T)+C4T^{C5}$$

C1:73.649、C2:-7,258.2、C3:−7.3037、C4:4.1653E−06、C5：2

ただし、圧力の単位はPa、適用温度範囲は、273.16～647.096[K]となっています。

上述のC1～C5パラメーターを各行毎に区切ってテキストに入力しておけば以下のような形で読み込むことができます。

①f=open('VPwat.txt')で同階層にあるVPwat.txtを開き、

②lines2=f.readlines()でテキストファイルを１行ずつlines2に格納していきます。

※readlines()以外にも、read()、readline()があります。使い分けは以下のリンクが分かりやすいので、ご参照ください。

read()

①f.close()でテキストファイルを閉じます。

④他の蒸気圧推算法でも対応できるように、VPの係数はMax7個まで格納できるようにしています。この際、Cに格納するデータが空だとエラーが出るため、空の部分には0を代入しています。

⑤C1～C7までデータを格納していきます。この際、リストからそのままデータを引っ張ってくると改行コードが邪魔であるのと、数値として認識されなかったため、.strip()で改行コードを除き、float()でそれぞれを数値化しています。

このように、係数は分離して別置きしておけば、推算式が変わっても呼び出し元の関数を使い分けるだけで対応することができます。

Excelファイルの読み取り

まずxlrdを読み込みます。

import xlrd

例えば、物性データを「PhysicalData.xlsx」として保存し、ここから必要な情報を読み取る場合、以下のように設定します。

Excelデータ

Pythonコード

簡単に各行のコードをご紹介します。

Excelファイルの読み込み

book = xlrd.open_workbook('PhysicalData.xlsx')

シートの読み込み（index(0)がsheet1）

sh1 = book.sheet_by_index(0)

セルの値の読み込み（例えば、sheet1のB2セルの値の読み取り）

sh1.cell(2, 2).value

純成分のzの求め方

3次状態方程式の解法1～２分法

3次状態方程式の解法２～はさみうち法

3次状態方程式の解法３～Cardano法とViete法

こちらのリンクを参照にしました。

上記のリンクによると、odeinitを使うと、硬い問題かどうかによって、採用手法を自動で切り替えてくれるとのこと。

odeを使うと、解法を自由に選択できるようなので、まずはodeを試してみます。

主要な解法を以下に記載します。

BDFの呼出し：ode(function).set_integrator('vode', method='bdf')

adamsの呼出し：ode(function).set_integrator('vode', method='adams')

lsodaの呼出し：ode(function).set_integrator('lsoda')

ode45の呼出し：ode(function).set_integrator('dopri5')

(functionの部分に関数定義します。）

例題として、stiffな問題の方が差異が確認しやすいので、Van Der Pol方程式で各手法を比較してみます。

比較前提

Time step:0.1[sec]、積分時間：2000[sec]、μ:1000

van der pol方程式は以下の通り。

$$\frac{d^2x}{dt^2}-\mu \left(1-x^2\right)\frac{dx}{dt}+x=0$$

Adamsは今回の条件では解くことがでなかったため、その他の計算手法の計算速度を比較してみると以下のようになります。

BDF：19.71[sec]、lsoda：0.14[sec]、ODE45：13.73

続いてodeintを試してみます。

odeint(function, X0, t)の形で定義します。（X0が初期値、tが時間定義になります）

計算時間わずか0.0164[sec]。

基本的にはまずodeintを使うという戦略は妥当そうです。

なお、上記VDPの解をmatplotlibで作図すると以下のようになります（BDF）。

モジュールを使わずに陽解法 モジュールを使わずにRunge-Kutta法

Python-Control

Pythonには制御検討用のモジュールが用意されていますので、こちらを使うと簡単に制御モデルの表現が可能です。

こちらでマニュアルを確認することができます。

control.matlabを使うとMatlab風に記載することができます。

　プロセス伝達関数 　

ステップ応答を見たい場合は、step()を使うことができます。

　

例えば、確認したいプロセスの伝達関数が1/(s^2+7s+10)で表せる場合、以下のように事前に伝達関数を定義し、それをstep関数に入れることで表現できます。

num = [0, 0, 1]
den = [1, 7, 10]
s = tf(num, den)
y, t = step(s, T=arange(0, 3, 0.1))

Bode線図の作成

Python-Controlを使うとBode線図の作成も容易にできます。

前節でご紹介した伝達関数sに対してBode線図を作成する場合は、以下のように１行のコードで済みます。

bode(s)

Nyquist線図の作成

Python-Controlを使うとBode線図同様、Nyquist線図も１行のコードで作成できます。

nyquist(s)

ただし、上記の伝達関数は1/(s^2+s+1)に対して書いたものです。

単一ページのHello Worldから、簡単な入力フォームを含むウェブアプリの作成方法をご紹介します。

単一ページのHello Worldから、csvファイルのアップロード、グラフ作成方法等ご紹介します。

Streamlitを使ったWebアプリの作成方法をご紹介します。

こちらのページに、簡単な説明を追加しました。

Pythonの学習に役立った講座、本をまとめていきます。

Udemyで参考になった講座

Udemyでは、ダントツで酒井さんの講座がお勧めです。

Courseraで参考になった講座

データの可視化にはData Visualization wiht Pythonがおすすめです。

Applied Machine Learning in Pythonがお勧めです。最終課題がちょっと大変でしたが、、

息子用に遊びで作ったコードです。意味は無いですが、対話型コードの参考として見てみてください。

While，If文の用法も確認でき、意外とこうした遊び感覚で作るプログラムの方が勉強になったりします。

※1 Shellでの入力要求はinput()の形で書きますが、そのままでは文字列として認識されてしまいます。

このため、数値として使用するためには、int(input())として一旦数値に変換する必要があります。

逆にprintで出力するには、文字列に変換する必要があるため、str()を使用します。