MIKU

7/11の活動について。

メトロノームの運動のモデルをやった。
式としては以下の通りになる。
ml\ddot{\theta}=-mg\sin\theta-\mu l\dot{\theta}+r
rは以下のように定める。
\omega>0,\theta_{1}<\theta<\theta_{2}のとき、r=r_{0}(1-\cos(\frac{2\pi(\theta-\theta_{1})}{\theta_{2}-\theta_{1}}))
\omega<0,-\theta_{2}<\theta<-\theta_{1}のとき、r=-r_{0}(1-\cos(\frac{2\pi(\theta-\theta_{1})}{\theta_{2}-\theta_{1}}))
それ以外のとき、r=0

以前2階の微分方程式を解くパッケージが見つからなかったが、
2階の微分方程式は1階連立微分方程式で表現できるので、odesolveで解けるみたいだった。

#パッケージを読み込む
library(odesolve)
#パラメーターの設定
g<-9.8
l<-1
u<-1
m<-1
parms<-c(g,l,u,m)
t<-c(0,0.01*(1:1000))
#加える外力
a<-0.25
b<-2*pi/0.09
c<-0.01
myf<-function(y){
if((y[1]>0.01)&(y[1]<0.1)&(y[2]>0)){myf<- a*(1-cos(b*(y[1]-c)))}else{
if((-0.01>y[1])&(-0.1<y[1])&(y[2]<0)){myf<- -a*(1-cos(b*(y[1]-c)))}else{myf<-0}}
return(myf)}
#微分方程式
dydt<-function(t,y,parms){
with(as.list(parms),{
dy1<- y[2]
dy2<- -g*sin(y[1])/l-u*y[2]/m+myf(y)
list(c(dy1,dy2))
})
}
#解いてみる
Y<-lsoda(c(0.09,0),t,dydt,parms)
plot(Y[,2],Y[,3],xlim=c(-0.1,0.1),ylim=c(-0.3,0.3))

図は以下の通りになった。