状態方程式へのfit をテンプレートにして作成

開始行:
**体積と圧力のデータをBirch&#8211;Murnaghan状態方程式（EO... 　３次Birch&#8211;Murnaghanの状態方程式(EOS)は、固体の状... 　nls()の()部分を変えれば、他の方程式のfitにも使える。 　以下に３次BM式にfitする例を示す。最初に体積データをx配... > x <- c(76.925555,76.3743575,75.82171,75.297905,74.784... 73.3523075,72.909195,72.47972,72.05904,71.65546,71.26874... 70.5173475,70.15515,69.7955775,69.45207,69.117005,68.786... 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... 　なお、多数のデータを打ち込むのは大変かつ間違いが入りや... > x = scan() > dat <- data.frame(x,y)　 > result <- nls(y ~ (3/2)*k0*((v0/x)^(7/3)-(v0/x)^(5/3))... ((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)) *... (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 ... 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をひっくり返している）... > 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,... > x.ax <- seq(65,78,0.01)　　# 以前この行が抜けてました... > points((3/2)*132.5*((76.93/x.ax)^(7/3)-(76.93/x.ax)^(5... ((76.93/x.ax)^(2/3)-1.0)),x.ax,col="blue",pch="-") > write.csv(dat,"Desktop/eos1.csv")　　# デスクトップへ... > q()　　# R終了

終了行:
**体積と圧力のデータをBirch&#8211;Murnaghan状態方程式（EO... 　３次Birch&#8211;Murnaghanの状態方程式(EOS)は、固体の状... 　nls()の()部分を変えれば、他の方程式のfitにも使える。 　以下に３次BM式にfitする例を示す。最初に体積データをx配... > x <- c(76.925555,76.3743575,75.82171,75.297905,74.784... 73.3523075,72.909195,72.47972,72.05904,71.65546,71.26874... 70.5173475,70.15515,69.7955775,69.45207,69.117005,68.786... 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... 　なお、多数のデータを打ち込むのは大変かつ間違いが入りや... > x = scan() > dat <- data.frame(x,y)　 > result <- nls(y ~ (3/2)*k0*((v0/x)^(7/3)-(v0/x)^(5/3))... ((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)) *... (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 ... 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をひっくり返している）... > 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,... > x.ax <- seq(65,78,0.01)　　# 以前この行が抜けてました... > points((3/2)*132.5*((76.93/x.ax)^(7/3)-(76.93/x.ax)^(5... ((76.93/x.ax)^(2/3)-1.0)),x.ax,col="blue",pch="-") > write.csv(dat,"Desktop/eos1.csv")　　# デスクトップへ... > q()　　# R終了

ページ名: