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

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

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

③グラフ描画も比較的容易

④外部ファイルの入出力が行いやすい

⑤変数宣言しないでガンガン書ける

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

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

・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

scipy.odeint

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法