計算統計力学 森北出版株式会社 Wm.G.Hoover

#p.9　ばねの運動方程式 d^2x/dt^2=-k*x/m
#2次差分を用いた場合(精度が悪い)

x=1 #初期置
xvec=c(x) #閾値ベクトル
tvec=c(0) #時間軸
ite=30 #計算回数
dt=1.2　#時間幅
k=5 #ばね定数
m=10 #質量

for(l in 1:ite){
if(l==1){
x=2*x-k*(dt^2)*x/m  
}else{
x=2*xvec[l]-k*(dt^2)*xvec[l]/m-xvec[l-1]  
}  
xvec=c(xvec,x)
tvec=c(tvec,l*dt)
}


plot(tvec,xvec,xlim=c(min(tvec),max(tvec)),ylim=c(min(xvec),max(xvec)),xlab="t",ylab="",col=2,type="o")

par(new=T)

plot(tvec,cos(sqrt(k/m)*tvec),xlim=c(min(tvec),max(tvec)),ylim=c(min(xvec),max(xvec)),xlab="t",ylab="",col=3,type="o")

#p.9 ばねの運動方程式 d^2x/dt^2=-k*x/m
#オイラー法

x=1 #初期置
y=0 #xの勾配の初期値
xvec=c(x) #閾値ベクトル
yvec=c(y)
tvec=c(0) #時間軸
ite=50 #計算回数
dt=0.1　#時間幅
k=5 #ばね定数
m=10 #質量

for(l in 1:ite){
x=x+dt*y
y=y-dt*k*x/m
xvec=c(xvec,x)
yvec=c(yvec,y)
tvec=c(tvec,l*dt)
}


plot(tvec,xvec,xlim=c(min(tvec),max(tvec)),ylim=c(min(xvec),max(xvec)),xlab="t",ylab="",col=2,type="o")

par(new=T)

plot(tvec,cos(sqrt(k/m)*tvec),xlim=c(min(tvec),max(tvec)),ylim=c(min(xvec),max(xvec)),xlab="t",ylab="",col=3,type="o")

#p.16 最小作用の原理
#円上を動く座標でラグランジアンを最小にするx,yの軌跡を求める
#オイラー法

x=1 #初期値x
y=0 #初期値y
X=0 #初期値xの勾配
Y=1 #初期値yの勾配
t=0 #初期時刻
lam=1 #初期値
r=2 #動く円上の半径

xvec=c(x) #閾値ベクトル
yvec=c(y) #閾値ベクトル
tvec=c(t) #時間軸

ite=1000
m=10 #質量
dt=0.01
Lvec=c() #ラグランジアン

for(l in 1:ite){
x=x+dt*X
y=y+dt*Y
x=r*round(x/sqrt(x^2+y^2),digit=4) #半径rの円上の運動に限定
y=r*round(y/sqrt(x^2+y^2),digit=4) #半径rの円上の運動に限定

X_pre=X
Y_pre=Y
X=X+lam*dt*x/m
Y=Y+lam*dt*y/m
dX=(X-X_pre)/dt
dY=(Y-Y_pre)/dt
lam=m*(dX+dY)/(x+y)
t=t+dt
L=m*(dX^2+dY^2)/2+lam*(x^2+y^2-r^2)/2

xvec=c(xvec,x)
yvec=c(yvec,y)
tvec=c(tvec,t)
Lvec=c(Lvec,L)
}

#ラグランジュ乗数の確認
m*dX/x
m*dY/y


plot(xvec,yvec)