MIKU
7/11の活動について。
メトロノームの運動のモデルをやった。
式としては以下の通りになる。
は以下のように定める。
のとき、
のとき、
それ以外のとき、
以前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))