体積と圧力のデータをBirch–Murnaghan状態方程式(EOS)にフィットする。

  • 3次Birch–Murnaghanの状態方程式(EOS)は、固体の状態方程式の1つで、地球物理分野ではよく使われる。ちょっとfitさせる必要があり、Rを試しに使ってみた。結構うまくいったので、メモしておく。 (訂正&追加:2015 11/12)
  • nls()の()部分を変えれば、他の方程式のfitにも使える。
  • 以下に3次BM式にfitする例を示す。最初に体積データをx配列に入れて、対応する圧力データをy配列に入れている。Rの非線形最小自乗法nlsを使ってV0, K0, K'を求めている。それらについては初期推定値が必要。V0は体積データから適当に、K'は初期値として4.0を使えばよい。K0は30 ~ 300 (GPa)くらいの値を取り、柔らかいものは小さく、硬いものは大きくする。変な値を入れないかぎり、収束する。
    > x <-  c(76.925555,76.3743575,75.82171,75.297905,74.784705,74.28808,73.8135775,
    73.3523075,72.909195,72.47972,72.05904,71.65546,71.26874,70.8880375,
    70.5173475,70.15515,69.7955775,69.45207,69.117005,68.7869625,68.46607,
    68.1508475,67.84432,67.5450325,67.248065,66.95757)
    > y <-  
    c(0.0001,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25)
  • なお、多数のデータを打ち込むのは大変かつ間違いが入りやすいので、数字をExcel等のテーブルからコピー&ペーストすることができる。そのためには、scan()を使う。使い方は、次行のように打ち込むと[1]と出るので、そこでデータをペーストする。データ終わりでリターンを押す。そのデータを直す時はfix(x)を使う。データがないところはNAを入れる。
    > x = scan() 
> dat <- data.frame(x,y) 
> result <- nls(y ~ (3/2)*k0*((v0/x)^(7/3)-(v0/x)^(5/3))*(1.0+(3/4)*(k1-4.0)*
((v0/x)^(2/3)-1.0)),dat,start=list(v0=77,k0=150,k1=4.0))
> summary(result) # 結果の表示

Formula: y ~ (3/2) * k0 * ((v0/x)^(7/3) - (v0/x)^(5/3)) * (1 + (3/4) * 
   (k1 - 4) * ((v0/x)^(2/3) - 1))

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
v0 7.693e+01  4.157e-03 18504.6   <2e-16 ***
k0 1.325e+02  2.243e-01   590.8   <2e-16 ***
k1 4.406e+00  2.029e-02   217.2   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.01265 on 23 degrees of freedom 

Number of iterations to convergence: 3 
Achieved convergence tolerance: 1.912e-08
  • 上記より正確な係数を知りたい時:
    > coefficients(result)

データとフィットした式のプロット

  • データ自体のプロット(都合上、xyをひっくり返している)。MacのRでの場合でしかチェックしてない。title()は軸タイトルにsuperscriptなど必要なので使っている。そんな必要なければ、xlab, ylabで直接指定すればよい。
    > plot(y,x,col="red",pch=19,xlab="",ylab="")
    > title(ylab=expression(paste("Volume (",&#197;^3,")")),line=2.5)
    > title(xlab=expression(paste("Pressure (GPa)")),line=2.5)
  • フィットした曲線のプロット:以下の式には上記で求めたV0,K0,K'を代入する。Vの範囲を最初にx.axに入れる。上記のデータをプロットしたウィンドウに上書きされるので、先に閉じないこと。Macの場合、このウィンドウをセーブするとpdfで保存される。また、以下はデータに応じて書き換えないといけない。
    > x.ax <- seq(65,78,0.01)  # 以前この行が抜けてました。
    > points((3/2)*132.5*((76.93/x.ax)^(7/3)-(76.93/x.ax)^(5/3))*(1.0+(3/4)*(0.406)*
    ((76.93/x.ax)^(2/3)-1.0)),x.ax,col="blue",pch="-")
    > write.csv(dat,"Desktop/eos1.csv")  # デスクトップへデータのcsvでの書き出し
    > q()  # R終了

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2018-01-15 (月) 10:56:15 (463d)