状態方程式へのfit の変更点

#author("2024-03-10T17:10:15+09:00","default:masami","masami")
RIGHT:&ref(https://www.misasa.okayama-u.ac.jp/~masami/images/OUOperatingGrant.png);
#author("2024-03-10T17:11:20+09:00;2024-03-10T17:10:15+09:00","default:masami","masami")
**体積と圧力のデータをBirch–Murnaghan状態方程式（EOS)にフィットする。 [#j6568f73]
-３次Birch–Murnaghanの状態方程式(EOS)は、固体の状態方程式の１つで、地球物理分野ではよく使われる。ちょっとfitさせる必要があり、Rを試しに使ってみた。結構うまくいったので、メモしておく。 (訂正＆追加：2015 11/12)
-nls()の()部分を変えれば、他の方程式のfitにも使える。
-以下に３次BM式にfitする例を示す。最初に体積データをx配列に入れて、対応する圧力データをy配列に入れている。Rの非線形最小自乗法nlsを使ってV0, K0, K'を求めている。それらについては初期推定値が必要。V0は体積データから適当に、K'は初期値として4.0を使えばよい。K0は30 ~ 300 (GPa)くらいの値を取り、柔らかいものは小さく、硬いものは大きくする。変な値を入れないかぎり、収束する。
　３次Birch–Murnaghanの状態方程式(EOS)は、固体の状態方程式の１つで、地球物理分野ではよく使われる。ちょっとfitさせる必要があり、Rを試しに使ってみた。結構うまくいったので、メモしておく。 (訂正＆追加：2015 11/12)
　nls()の()部分を変えれば、他の方程式のfitにも使える。
　以下に３次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を入れる。
　なお、多数のデータを打ち込むのは大変かつ間違いが入りやすいので、数字を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)
**データとフィットした式のプロット [#w0b67236]
-データ自体のプロット（都合上、xyをひっくり返している）。MacのRでの場合でしかチェックしてない。title()は軸タイトルにsuperscriptなど必要なので使っている。そんな必要なければ、xlab, ylabで直接指定すればよい。２番目のところはAngstrom記号を使ってますが、変なコードで表示される時があります。その場合は適当に変えてください。
　データ自体のプロット（都合上、xyをひっくり返している）。MacのRでの場合でしかチェックしてない。title()は軸タイトルにsuperscriptなど必要なので使っている。そんな必要なければ、xlab, ylabで直接指定すればよい。２番目のところはAngstrom記号を使ってますが、変なコードで表示される時があります。その場合は適当に変えてください。
 > plot(y,x,col="red",pch=19,xlab="",ylab="")
 > title(ylab=expression(paste("Volume (Å^3,")")),line=2.5)
 > title(xlab=expression(paste("Pressure (GPa)")),line=2.5)
-フィットした曲線のプロット：以下の式には上記で求めたV0,K0,K'を代入する。Vの範囲を最初にx.axに入れる。上記のデータをプロットしたウィンドウに上書きされるので、先に閉じないこと。Macの場合、このウィンドウをセーブするとpdfで保存される。また、以下はデータに応じて書き換えないといけない。
　フィットした曲線のプロット：以下の式には上記で求めた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)*(4.406-4.0)*
 ((76.93/x.ax)^(2/3)-1.0)),x.ax,col="blue",pch="-")
 > write.csv(dat,"Desktop/eos1.csv")　　# デスクトップへデータのcsvでの書き出し
 > q()　　# R終了