Python

化学工学、数値計算に役立つTipsを紹介します。

まだコンテンツは少ないですが、今後状態方程式法の解法、常微分方程式の解法等、数値計算、制御のTipsを追加していく予定です。

事前準備

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

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

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

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

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

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

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

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

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

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

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

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

よく使うモジュール(Numpy, Scipy, 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で計算する場合、非常に短いコードで済みます。

コード例            出力例

SimpleMatrix

このように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、活量係数をまとめてプロットしたいという場合、以下のようにまとめて表示させることが可能です。

NRTLplot

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

subplots 3dplot出力

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

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

MatplotlibGraph 三角座標のプロット

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

活用例

例えば、以下の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パラメーターを各行毎に区切ってテキストに入力しておけば以下のような形で読み込むことができます。

VPreadCode

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

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

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

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

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

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

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

    Excelファイルの読み取り

    まずxlrdを読み込みます。

    import xlrd

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

    Excelデータ

    ExcelData

    Pythonコード

    xlrdread

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

    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~2分法 3次状態方程式の解法2~はさみうち法 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は今回の条件では解くことがでなかったため、その他の計算手法の計算速度を比較してみると以下のようになります。

    ODE計算時間比較

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

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

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

    odeinit

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

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

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

    VDPplot

    準備中

    準備中

       

    制御の基本

     プロセス伝達関数
    Bode線図の作成 Nyquist線図の作成

    おまけ

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

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

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

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

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

    当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。

    Page Top