063_Vibration2.mw

機械振動論 

2質点系 

2質点系 

次に2質点系のモデルを考えてみましょう。 

 

1質点系の一輪車のモデルの上にもうひとつバネと物体をつないだものです。 

 

Mapleを初期化します。 

> restart:
 

プロットのパッケージも用意しておきます。 

> with(plots):with(plottools):
 

Warning, the name changecoords has been redefined 

Warning, the assigned name arrow now has a global binding 

運動方程式は二つの方程式系(連立微分方程式)となります。 

> eq1:=M1*(diff(x1(t),t$2)+G)+c1*diff(x1(t),t)+k1*(x1(t)-g(t))=0;
 

M1*((diff(diff(x1(t), t), t))+G)+c1*(diff(x1(t), t))+k1*(x1(t)-g(t)) = 0
M1*((diff(diff(x1(t), t), t))+G)+c1*(diff(x1(t), t))+k1*(x1(t)-g(t)) = 0
 

> eq2:=M2*(diff(x2(t),t$2)+G)+c2*diff(x2(t),t)+k2*(x2(t)-x1(t)-l2)=0;
 

M2*((diff(diff(x2(t), t), t))+G)+c2*(diff(x2(t), t))+k2*(x2(t)-x1(t)-l2) = 0
M2*((diff(diff(x2(t), t), t))+G)+c2*(diff(x2(t), t))+k2*(x2(t)-x1(t)-l2) = 0
 

方程式をeq1とeq2に保存します。 

 

質量はそれぞれM1、M2とします。 

x1(t)はM1の高さ、x2(t)はM2の時間tでの高さです。 

減衰定数はc1,c2とします。 

バネ定数はc3とします。 

上のバネの長さ(質点間の距離)l2とします。 

 

これもMapleのdsolveにより数式的に解くことができます。 

> dsolve({eq1,eq2},{x1(t),x2(t)});
 

{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
{x2(t) = (_C2*exp((-1/2*c2+1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+_C1*exp((-1/2*c2-1/2*(c2^2-4*k2*M2)^(1/2))*t/M2)*(c2^2-4*k2*M2)^(1/2)+Int(exp(-1/2*t*(-M1*c2+M1*(c2^2-4*k2*M2)^(1/2)+M2*...
 

このように線形方程式なのでMapleのdsolveは解を見つけてくれますが、これではあまりにも複雑すぎて意味がありません。 

 

具体的な値を設定してもっと簡単にしてみましょう。 

 

質量を定義します。 

> M1:=1;M2:=1.2;
 

1 

1.2 

重力加速度はMKSでは9.8です。 

> G:=9.8;
 

9.8 

減衰定数c1、c2は0.1と0.2にします。 

> c1:=0.1;c2:=0.2;
 

.1 

.2 

バネ定数k1、k2は100と200とします。 

> k1:=100;k2:=200;
 

100 

200 

質点間の距離l2は2とします。 

> l2:=2;
 

2 

地面の変化(強制振動)は単純な正弦関数で次のようにします。 

> g:=t->sin(t)/5;
 

proc (t) options operator, arrow; 1/5*sin(t) end proc 

これで運動方程式は次のようになります。 

> eq1;eq2;
 

(diff(diff(x1(t), t), t))+9.8+.1*(diff(x1(t), t))+100*x1(t)-20*sin(t) = 0 

1.2*(diff(diff(x2(t), t), t))-388.24+.2*(diff(x2(t), t))+200*x2(t)-200*x1(t) = 0
1.2*(diff(diff(x2(t), t), t))-388.24+.2*(diff(x2(t), t))+200*x2(t)-200*x1(t) = 0
 

dsolveで解きます。M2の初期値の位置(高さ)はM1の上l2のところです。その他の初期値は0とします。 

> ans:=dsolve({eq1,eq2,x1(0)=0,x2(0)=l2,D(x1)(0)=0,D(x2)(0)=0},{x1(t),x2(t)});
 

{x2(t) = 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*...
{x2(t) = 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*...
{x2(t) = 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*...
{x2(t) = 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*...
{x2(t) = 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*...
{x2(t) = 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*...
{x2(t) = 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*...
{x2(t) = 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*...
 

x1(t)の解です。 

> ans1:=rhs(ans[1]);
 

1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*39999^(1/...
1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*39999^(1/...
1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*39999^(1/...
1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*39999^(1/...
1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/78406119798000000*exp(-1/20*t)*sin(1/20*39999^(1/...
 

x2(t)の解です。 

> ans2:=rhs(ans[2]);
 

-1931875051/19601529949500*exp(-1/20*t)*sin(1/20*39999^(1/2)*t)*39999^(1/2)+48124949/490050500*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-49/500-200/980101*cos(t)+198000/980101*sin(t)
-1931875051/19601529949500*exp(-1/20*t)*sin(1/20*39999^(1/2)*t)*39999^(1/2)+48124949/490050500*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-49/500-200/980101*cos(t)+198000/980101*sin(t)
-1931875051/19601529949500*exp(-1/20*t)*sin(1/20*39999^(1/2)*t)*39999^(1/2)+48124949/490050500*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-49/500-200/980101*cos(t)+198000/980101*sin(t)
 

unapplyでx1とx2を関数化します。 

> x1:=unapply(ans1,t);x2:=unapply(ans2,t);
 

proc (t) options operator, arrow; 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/784061197980000...
proc (t) options operator, arrow; 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/784061197980000...
proc (t) options operator, arrow; 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/784061197980000...
proc (t) options operator, arrow; 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/784061197980000...
proc (t) options operator, arrow; 1152/625-396800000/968376051737*cos(t)+196811800000/968376051737*sin(t)+482191365051/1960202000000*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-18355309654949/784061197980000...
 

proc (t) options operator, arrow; -1931875051/19601529949500*exp(-1/20*t)*sin(1/20*39999^(1/2)*t)*39999^(1/2)+48124949/490050500*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-49/500-200/980101*cos(t)+198000/98...
proc (t) options operator, arrow; -1931875051/19601529949500*exp(-1/20*t)*sin(1/20*39999^(1/2)*t)*39999^(1/2)+48124949/490050500*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-49/500-200/980101*cos(t)+198000/98...
proc (t) options operator, arrow; -1931875051/19601529949500*exp(-1/20*t)*sin(1/20*39999^(1/2)*t)*39999^(1/2)+48124949/490050500*exp(-1/20*t)*cos(1/20*39999^(1/2)*t)-49/500-200/980101*cos(t)+198000/98...
 

2つの高さの関数(解)をプロットしてみましょう。 

> plot({x1(t),x2(t)},t=0..4);
 

Plot 

下(緑)がM1、上(赤)の方がM2の高さです。 

 

これをまたアニメーションにしてみましょう。 

> p:=[]:
 

pはアニメーションのフレームを追加していくための変数で、最初は空にしておきます。 

> r:=1:
 

rは車輪の半径です。 

> c:=cuboid([-0.2,-0.2,-0.2],[0.2,0.2,1.2],color=green):
 

cはタイヤにつながっている部分を直方体(cuboid)で作成しておきます。 

対角の2点の座標で定義します。色は緑(green)にします。 

これをtranslateでシフトしてアニメーションします。 

> M1:=cuboid([-0.5,-0.5,0],[0.5,0.5,1],color=red):
 

> M2:=cuboid([-0.5,-0.5,0],[0.5,0.5,1],color=green):
 

M1の色は赤、M2の色は緑にします。 

 

以下のforループで時間(tが0から1秒まで)ごとにアニメーションのフレームを作成します。 

> for t from 0 by 0.05 to 5 do
 h:=g(t):
 pp:=plot3d([x,y,g(x+t)],x=-Pi/2..Pi/2,y=-1..1,grid=[15,12],color=gray):
 pp:=pp,translate(rotate(cylinder([0,0,0],r,0.2,color=blue),Pi/2,t*2*Pi,0),0,-0.1,h+r):
 pp:=pp,translate(c,0,0,h+r):
 H:=x1(t):
 L:=H-h+1:
 tt:=tubeplot([cos(u*8*Pi)/4,sin(u*8*Pi)/4,u*L],u=0..1,radius=0.05,color=yellow):
 pp:=pp,translate(tt,0,0,h+2*r+0.2):
 pp:=pp,translate(M1,0,0,h+2*r+0.2+L):
 tt:=tubeplot([cos(u*8*Pi)/4,sin(u*8*Pi)/4,u*(x2(t)-x1(t))],u=0..1,radius=0.05,color=cyan):
 pp:=pp,translate(tt,0,0,h+2*r+0.2+L+0.5):
 pp:=pp,translate(M2,0,0,h+2*r+0.2+L+x2(t)-x1(t)):
 p:=p,display(pp):
end:
 

できたフレームをアニメーション表示します。 

> display(p,insequence=true,scaling=constrained);
 

Plot 

>
 

行列表現 

この2質点系の運動方程式は行列で表現することができます。 

 

またMapleを初期化します。 

> restart:
 

列(column)ベクトルを表現するために線形代数パッケージを呼び出します。 

> with(LinearAlgebra):
 

二階微分のベクトルddxを用意します。 

> ddx:=Vector([diff(x1(t),t$2),diff(x2(t),t$2)]);
 

Vector[column](%id = 677380016) 

一階微分のベクトルdxを用意します。 

> dx:=Vector([diff(x1(t),t),diff(x2(t),t)]);
 

Vector[column](%id = 677442664) 

位置(高さ)のベクトルを用意します。 

> xx:=Vector([x1(t),x2(t)]);
 

Vector[column](%id = 677442784) 

右辺の入力のベクトルffを定義します。 

> ff:=Vector([k1*g(t),k2*l2]);
 

Vector[column](%id = 677442824) 

これで2質点系の運動方程式は次のようになります。ここで&*はMapleでは行列の積を表し 

通常のスカラーの積*と区別します(非可換積)。 

行列の積では A &* B ≠ B &* A です。 

> <<M1|0>,<0|M2>>&*ddx+<<c1|0>,<0|c2>>&*dx+<<k1|0>,<-k2|k2>>&*xx=ff;
 

`&*`(Matrix(%id = 677724180), Vector[column](%id = 677380016))+`&*`(Matrix(%id = 677765980), Vector[column](%id = 677442664))+`&*`(Matrix(%id = 677804724), Vector[column](%id = 677442784)) = Vector[co...
`&*`(Matrix(%id = 677724180), Vector[column](%id = 677380016))+`&*`(Matrix(%id = 677765980), Vector[column](%id = 677442664))+`&*`(Matrix(%id = 677804724), Vector[column](%id = 677442784)) = Vector[co...
 

それぞれ係数の行列をM,C,Kで表し、ベクトル関数(x1(t),x2(t))をx(t)で表し 

右辺の入力ベクトル関数をf(t)とすると次のようにな一般化された行列表現なります。 

> M*diff(x(t),t$2)+C*diff(x(t),t)+K*x(t)=f(t);
 

M*(diff(diff(x(t), t), t))+C*(diff(x(t), t))+K*x(t) = f(t) 

>
 

行列表現する意味は多質点系のシステムになったときにも一般化できるからです。