067_Vibration6.mw

機械振動論 

平面振動 

質点を平面上にならべて平面状の系の振動をシミュレーションします。 

平面上に並べた質点系 

まずは解りやすいように3D表示しましょう。 

> restart:
 

> with(plots):with(plottools):
 

Warning, the name changecoords has been redefined 

Warning, the assigned name arrow now has a global binding 

質点を表す立方体を用意します。 

> cb:=cuboid([0,0,0],[1,1,1]):
 

表示してみます。 

> display(cb,scaling=constrained);
 

Plot 

スプリングを用意します。 

> spr:=tubeplot([cos(4*Pi*t)/5,sin(4*Pi*t)/5,t],t=0..2,radius=0.1):
 

回転させておきます。 

> spr:=rotate(spr,Pi/2,0,0):
 

表示してみます。 

> display(spr,scaling=constrained);
 

Plot 

> n:=5;
 

5 

n^2 の質点の配置です。 

オブジェクトを保存する変数を空にします。 

> p:=[];
 

[] 

質点とスプリングを配置します。 

> for i to n do
 for j to n do
   p:=p,translate(cb,3*(i-1),3*(j-1),0):
   p:=p,translate(spr,3*(i-1)+0.5,3*(j-1)+1,0.5):
   p:=p,translate(rotate(spr,0,0,Pi/2),3*(i-1)+1,3*(j-1)+0.5,0.5):
 end:
end:
 

端の質点とスプリングを追加します。 

> for i to n do
 p:=p,translate(cb,3*(i-1),-3,0):
 p:=p,translate(cb,-3,3*(i-1),0):
 p:=p,translate(cb,3*(i-1),3*n,0):
 p:=p,translate(cb,3*n,3*(i-1),0):
 p:=p,translate(spr,3*(i-1)+0.5,-2,0.5):
 p:=p,translate(rotate(spr,0,0,Pi/2),-2,3*(i-1)+0.5,0.5):
end:
 

表示します。 

> display(p,scaling=constrained);
 

Plot 

>
 

運動方程式 

この系の運動方程式を立てます。 

> n:=5;
 

5 

25 点です。
運動方程式
 

> for i to n do
 for j to n do
   eq[i,j]:=m*diff(x[i,j](t),t$2)+c*diff(x[i,j](t),t)+k*(4*x[i,j](t)-x[i+1,j](t)-x[i-1,j](t)-x[i,j+1](t)-x[i,j-1](t))=0;
 end;
end;
 

パラメータを決定します。 

> m:=1;c:=0.5;k:=100;
 

1 

.5 

100 

方程式のセットを構成します。evalfで数値化します。 

> equ:={}:
 

> for i to n do
 for j to n do
   equ:=equ union {evalf(eq[i,j])}:
 end:
end:
 

端の質点は固定します。すべて値0の関数です。 

> for i to n do
 x[i,0]:=t->0;
 x[i,n+1]:=t->0;
 x[0,i]:=t->0;
 x[n+1,i]:=t->0;
end:
 

初期値のセット(高さと初速)を空にします。 

> init:={};init1:={};
 

{} 

{} 

初期値(t=0)をすべて0にセットします。 

> for i to n do
 for j to n do
   init:=init union {x[i,j](0)=0}:
   init1:=init1 union {D(x[i,j])(0)=0}:
 end:
end:
 

(2,2)の質点の初速の初期値を外します。 

> init1:=init1 minus {D(x[2,2])(0)=0}:
 

あらためて(2,2)の初速を0.01に設定し追加します。 

> init1:=init1 union {D(x[2,2])(0)=0.01};
 

{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
{(D(x[2, 4]))(0) = 0, (D(x[5, 5]))(0) = 0, (D(x[2, 5]))(0) = 0, (D(x[1, 1]))(0) = 0, (D(x[1, 2]))(0) = 0, (D(x[3, 1]))(0) = 0, (D(x[1, 3]))(0) = 0, (D(x[3, 2]))(0) = 0, (D(x[5, 3]))(0) = 0, (D(x[1, 4]...
 

解の関数群(n^2)を決定します。 

> funcs:={}:
 

> for i to n do
 for j to n do
   funcs:=funcs union {x[i,j](t)}
 end:
end:
 

解のサンプルの個数を決めます。 

> m:=50;
 

50 

解のサンプルの時間刻みをきめます。0秒からπ秒まで 

> d:=Pi/m;
 

1/50*Pi 

サンプル時間のセット 

> samp:=array([seq(j*d,j=0..m)]):
 

dsolceで微分方程式を解きます。 

> ans:=dsolve(equ union init union init1,funcs,numeric,output=samp):
 

解の関数リストです。 

> seq(ans[1,1][j],j=1..2*n*n+1);
 

t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
t, x[1, 1](t), diff(x[1, 1](t), t), x[1, 2](t), diff(x[1, 2](t), t), x[1, 3](t), diff(x[1, 3](t), t), x[1, 4](t), diff(x[1, 4](t), t), x[1, 5](t), diff(x[1, 5](t), t), x[2, 1](t), diff(x[2, 1](t), t),...
 

(1,1)の解のプロット 

> listplot([seq(ans[2,1][j,2],j=1..m+1)]);
 

Plot 

(1,2)の解のプロット 

> listplot([seq(ans[2,1][j,4],j=1..m+1)]);
 

Plot 

>
 

アニメーション 

解を使ってアニメーションします。 

 

まず質点を表す平板を用意します。 

> cb:=cuboid([0,0,0],[1,1,0.1]):
 

フレームを保存する変数pを空にします。 

> p:=[];
 

[] 

時間をすすめながらフレームを追加していきます。 

変化は1000倍に誇張して表現しています。 

> for k from 0 to m do
 pp:=[]:
 for i to n do
   for j to n do
     h:=ans[2,1][k+1,2*((i-1)*n+j)]:
     pp:=pp,translate(cb,i,j,h*1000):
   end:
 end:
 p:=p,display(pp):
end:
 

アニメーションを表示します。 

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

Plot 

>
 

もう一つのケース 

(2,2)の初速の初期値を外します。 

> init1:=init1 minus {D(x[2,2])(0)=0.01}:
 

初速を0に戻します。 

> init1:=init1 union {D(x[2,2])(0)=0}:
 

(3,3)の初期値(初速)を外します。 

> init1:=init1 minus {D(x[3,3])(0)=0}:
 

初速を0.01にします。 

> init1:=init1 union {D(x[3,3])(0)=0.01}:
 

dsolveで解きます。 

> ans:=dsolve(equ union init union init1,funcs,numeric,output=samp):
 

アニメーションを作成します。pを空にします。 

> p:=[];
 

[] 

時間をすすめながらフレームを追加していきます。 

変化は1000倍に誇張して表現しています。 

> for k from 0 to m do
 pp:=[]:
 for i to n do
   for j to n do
     h:=ans[2,1][k+1,2*((i-1)*n+j)]:
     pp:=pp,translate(cb,i,j,h*1000):
   end:
 end:
 p:=p,display(pp):
end:
 

アニメーションを表示します。 

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

Plot