統計力学   裳華房　阿部豊先生

物理の勉強が必要かと。

```r

#小正準集団(ミクロカノニカル分布) p.27
#粒子数N、体積V、エネルギーEが一定

J=5
ep=sample(1:5,J,replace=T)#各1粒子のとりうるエネルギー
del=sample(1:5,J,replace=T)#とりうる状態の個数
nj=sample(3:7,J,replace=T)#全粒子数Nの割り当て
N=sum(nj)
E=sum(ep*nj)#全エネルギー
kB=1.381*10^(-23)


f=function(x){
val=cbind(rep(1,J),ep)%*%x
fval=c(sum(del*exp(-val))-N,sum(del*ep*exp(-val))-E)
return(fval)
}


#ニュートン法
X=rep(1,2)

ite=120000
eta=10^(-4)
h=10^(-5)

for(l in 1:ite){
dF=array(0,dim=c(length(X),length(X))) 
for(j in 1:length(X)){
Xsub=X
Xsub[j]=Xsub[j]+h
dF[,j]=(f(Xsub)-f(X))/(h)
}  
 
X=X-eta*solve(dF)%*%f(X)
#print(sum(f(X)))   
}  

#ラグランジュ定数
alpha=X[1]
beta=X[2]

#温度T(3.16)?
t=1/(kB*beta)

#W_nj(3.3)を最大にするnj
nj_star=del*exp(-cbind(rep(1,J),ep)%*%X)
E_star=sum(ep*nj_star)


```

```r

#正準集団(カノニカル分布) p.28
#温度Tが指定された扱い

J=7
Mj=sample(1:5,J,replace=T)
M=sum(Mj)
Ej=sample(1:5,J,replace=T)
E=sum(Mj*Ej)
kB=1.381*10^(1)
t=298.15#室温が25℃の場合

canonical_dis=function(x){
return(exp(-x*Ej)/sum(exp(-x*Ej)))  
}

logW=function(x){
Mj2=M*canonical_dis(x)
return(M*log(M)-sum(Mj2*log(Mj2)))
}

Energy=function(x){#入力変数は温度
return(sum(Ej*exp(-x*Ej))/sum(exp(-x*Ej)))
}

#ヘルムホルツの自由エネルギー
f=function(x){
return(-kB*t*log(sum(exp(-x*Ej))))
}


beta=1/(kB*t)
h=10^(-5)
p_Ej=canonical_dis(beta)
E_star=Energy(beta)
VE_star=sum(Ej*Ej*exp(-beta*Ej))/sum(exp(-beta*Ej))-(sum(Ej*exp(-beta*Ej))/sum(exp(-beta*Ej)))^2
Cv=-kB*(beta^2)*(Energy(beta+h)-Energy(beta))/(h)#比熱
Mj_star=M*p_Ej  


#エントロピー
S=kB*(beta^2)*(f(beta+h)-f(beta))/h

```

```r

#大正準集団(グランドカノニカル分布) p.34
#

I=5
J=6
Eij=matrix(sample(1:5,I*J,replace=T),nrow=I)
Ni=sample(1:5,I,replace=T)
Mij=matrix(sample(1:5,I*J,replace=T),nrow=I)
M=sum(Mij)
E0=sum(Eij*Mij)
N0=sum(Mij*Ni)

f=function(x){
Mij2=M*exp(-x[1]*Eij-x[2]*Ni)/sum(exp(-x[1]*Eij-x[2]*Ni))
fval=c(sum(Eij*Mij2)-E0,sum(Mij2*Ni)-N0)
return(fval)
}


#ニュートン法
X=rep(1,2)
ite=100
eta=10^(-1)
h=10^(-5)

for(l in 1:ite){

dF=array(0,dim=c(length(X),length(X))) 
for(j in 1:length(X)){
Xsub=X
Xsub[j]=Xsub[j]+h
dF[,j]=(f(Xsub)-f(X))/(h)
}  
 
X=X-eta*solve(dF)%*%f(X)
print(sum(f(X)))   
}  

#ラグランジュ未定乗数
beta=X[1]
r=X[2]

Zg=sum(exp(-beta*Eij-r*Ni))#規格化定数
Pij_star=exp(-beta*Eij-r*Ni)/Zg#logWを最大にするように計算された確率
Mij_star=M*Pij_star#logWを最大にするように計算されたコピーの数

#エネルギー平均
E_ave=sum(Eij*exp(-beta*Eij-r*Ni))/Zg

#粒子数の平均
N_ave=sum(exp(-beta*Eij-r*Ni)*Ni)/Zg

#付加条件の確認
sum(Eij*Mij_star)#E0に等しい
sum(Mij_star*Ni)#N0に等しい


```

```r

#p.39
#(4.4)の計算

h=6.626*10^(-1)#プランク定数。6.626*10^(-34)だと計算しにくい。
I=10#慣性モーメント。とりあえずおいとく
beta=10^(2)
hh=10^(-5)
sig=beta*h^2/(2*I)

if(beta*h^2/I<1){

fl=function(x){
El=x*(x+1)*h^2/(2*I)  
return((2*x+1)*exp(-beta*El))
}

#オイラーマクローリン総和公式
Z1=integrate(fl,0,Inf)$value+fl(0)/2-(fl(0+hh)-fl(0))/(12*hh)  

flh=function(x){
El=x*(x+1)*h^2/(2*I)  
return((2*x+1)*exp(-(beta+hh)*El))
}

Z1h=integrate(flh,0,Inf)$value+flh(0)/2-(flh(0+hh)-flh(0))/(12*hh)
E=-(log(Z1h)-log(Z1))/hh


}else{

Z1=function(b){
return(1+3*exp(-b*h^2/I))  
}

E=function(b){
  return(-(log(Z1(b+hh))-log(Z1(b)))/hh)
}

Cv=-kB*(beta^2)*(E(beta+hh)-E(beta))/hh

}

```

```r

#Bose統計とFermi統計 p.62

#Bose統計
J=5
ep=sample(1:10,J,replace=T)#各1粒子のとりうるエネルギー
del=sample(1:5,J,replace=T)#とりうる状態の個数
nj=sample(3:7,J,replace=T)#全粒子数Nの割り当て
N=sum(nj)
E=sum(ep*nj)#全エネルギー

f=function(x){
val=cbind(rep(1,J),ep)%*%x
fval=c(sum(ep*del/(exp(val)-1))-E,sum(del/(exp(val)-1))-N)
return(fval)
}


#ニュートン法
X=rep(0.1,2)
ite=2000
eta=10^(-2)
h=10^(-5)
hh=0.01

for(l in 1:ite){

dF=array(0,dim=c(length(X),length(X))) 
for(j in 1:length(X)){
Xsub=X
Xsub[j]=Xsub[j]+h
dF[,j]=(f(Xsub)-f(X))/(h)
}  
 
X=X-eta*solve(dF)%*%f(X)
print(sum(f(X)))   
}  

#ラグランジュ未定乗数
alpha=X[1]
beta=X[2]


nj_star=del/(exp(cbind(rep(1,J),ep)%*%X)-1)
E_star=sum(ep*del/(exp(cbind(rep(1,J),ep)%*%X)-1))




```

```r

#Bose統計とFermi統計 p.62

#Fermi統計
J=5
ep=sample(1:10,J,replace=T)#各1粒子のとりうるエネルギー
del=sample(1:10,J,replace=T)#とりうる状態の個数
nj=sample(0:1,J,replace=T)#全粒子数Nの割り当て
N=sum(nj)
E=sum(ep*nj)#全エネルギー

f=function(x){
#val=x[1]*(ep-x[2])
#fval=c(sum(ep*del/(exp(val)+1))-E,sum(del/(exp(val)+1))-N)
val=cbind(rep(1,J),ep)%*%x
fval=c(sum(ep*del/(exp(val)+1))-E,sum(del/(exp(val)+1))-N)
return(fval)
}


#ニュートン法
X=c(1,1)
ite=10^3
eta=10^(-2)
h=10^(-5)
hh=0.01

for(l in 1:ite){

dF=array(0,dim=c(length(X),length(X))) 
for(j in 1:length(X)){
Xsub=X
Xsub[j]=Xsub[j]+h
dF[,j]=(f(Xsub)-f(X))/(h)
}  
 
X=X-eta*solve(dF)%*%f(X)
print(sum(f(X)))   
}  

#ラグランジュ未定乗数
alpha=X[1]
beta=X[2]


nj_star=del/(exp(cbind(rep(1,J),ep)%*%X)+1)
E_star=sum(ep*del/(exp(cbind(rep(1,J),ep)%*%X)+1))



```