11
24

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

[Pythonによる科学・技術計算] 非線形関数によるフィッティング, 状態方程式, scipy

Last updated at Posted at 2017-07-18

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

t.png


##[付録] 使った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
11
24
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
11
24

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?