scipy.optimizeのcurvefitメソッドを用いて， XYデータを非線形関数でフィットし， 最適なフィッティングパラーメータを得る。
目的: XYデータとして， 物質の体積(V)と圧力(P)データを準備し，
それを固体の体積・圧力の関係を表す状態方程式の一つとして広く使われている3次のバーチ・マーナガン状態方程式へとフィットすることでフィッティングパラメータを求める。
この方程式は以下の形を持つ。
P(V)=\frac{3}{2}K_0[(\frac{V}{V_0})^{-7/3}-(\frac{V}{V_0})^{-5/3}][1+\frac{3}{4}(K_d-4)((\frac{V}{V_0})^{-2/3}-1)]
$K_0, V_0, K_d$がフィッティングパラメータである。
また，X軸がV, y軸がPである。
問題によって，コード中のフィット関数形の指定を変えれば，あらゆる関数形で最適なフィットパラメータを得ることができる。また，線形でのフィットももちろん可能。特殊な状況を除き，万能なコードとなる。
import numpy as np
import scipy.optimize
import matplotlib.pylab as plt
"""
非線形関数によるカーブ・フィッティング
例: 圧力 vs 体積の関係を3次のバーチ・マーナガン方程式でフィッティングする
"""
x_list=[] # x_listを定義 (空のリストを作成)
y_list=[] # y_listを定義 (空のリストを作成)
f=open('EoS.dat','rt') # プロットしたいデータが入っているファイルをr(読み込み) t(テキスト)モードで読み込む
## データを読み込み，x_listとy_listに値を格納する
for line in f:
data = line[:-1].split(' ')
x_list.append(float(data[0]))
y_list.append(float(data[1]))
##
plt.plot(x_list, y_list, 'o',label='Raw data')
#フィッティング関数の指定
def fitfunc(V, V0, K0, K0d ):
return (3.0/2.0)*K0*((V/V0)**(-7.0/3.0)-(V/V0)**(-5.0/3.0))*(1.0+(3.0/4.0)*(K0d-4.0)*((V/V0)**(-2.0/3.0)-1.0))
para_ini =[200.0, 100.0, 4.0] # フィッティングパラメータの初期条件
para_opt, cov = scipy.optimize.curve_fit(fitfunc, x_list, y_list, para_ini)#フィッティングパラメータの最適化
#得られたフィッティングパラメータの表示
print ("Fitted V0 =", str(para_opt[0]), "(Å^3)")
print ("Fitted K0 =", str(para_opt[1]), "(GPa)")
print ("Fitted K0' =", str(para_opt[2]))
##
# for plot
V0=para_opt[0]
K0=para_opt[1]
K0d=para_opt[2]
V1 = np.arange(min(x_list) - 1, max(x_list) + 1, 0.1)
y_fit = (3.0/2.0)*K0*((V1/V0)**(-7/3)-(V1/V0)**(-5.0/3.0))*(1.0+(3.0/4.0)*(K0d-4.0)*((V1/V0)**(-2.0/3.0)-1.0))
plt.plot(V1, y_fit, color='RED',linewidth=2.0, label='Fitted curve')
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.grid(True)
plt.title("Curve fitting for eqaution of states")
plt.legend(loc='upper right')
plt.xlabel('Volume (Å^3)', fontsize=18)
plt.ylabel('Pressure (GPa)', fontsize=18)
plt.show()
結果
Fitted V0 = 163.627068001 (Å^3)
Fitted K0 = 225.52189108 (GPa)
Fitted K0' = 4.06160382625
[付録] 使ったEoS.datファイルの中身
164.26 0.0
156.05 11.9
147.83 27.9
146.19 31.7
144.55 35.7
142.91 39.9
141.26 44.4
139.62 49.2
137.98 54.2
136.34 59.5
134.69 65.1
133.05 71.0
131.41 77.3
129.77 83.9
128.12 90.9
126.48 98.3
124.84 106.2
123.20 114.4
121.55 123.2
119.91 132.5
118.27 142.3
116.62 152.6
114.98 163.6
113.34 175.2
111.70 187.5
110.05 200.6