064_Vibration3.mw

機械振動論

多質点系

ここでは更に多質点のシミュレーションに挑戦します。
今度は質点を上に積み上げた建物の横揺れのシミュレーションを行います。

Mapleを初期化しておきます。
> restart
 

「お知らせ」
この事例には誤りがあったため修正しております。
ユーザ様のご指摘で減衰の項も階ごとの差にすべきところが差になっていませんでした。
他にも誤りがある可能性がございます。
お気づきの方はぜひご指摘頂きますようお願いいたします。

減衰項
> c[1]*((diff(x[2](t), t))-(diff(x[1](t), t)))

(Typesetting:-mprintslash)([c[1]*((diff(x[2](t), t))-(diff(x[1](t), t)))], [c[1]*((diff(x[2](t), t))-(diff(x[1](t), t)))]) (1)

3質点

2質点の時と同様に今度は三つの運動方程式となります。
3階の建物とし各階を1質点としてモデリングします。

一階の運動方程式です。今回は縦のバネではなく横のバネでつながっているとします。
上下の運動は無視します。

> eq[1]:=m[1]*diff(x[1](t),t$2)+c[1]*(diff(x[1](t),t)-diff(x[0](t),t))+k[1]*(x[1](t)-x[0](t))=0;

(Typesetting:-mprintslash)([eq[1] := m[1]*(diff(x[1](t), `$`(t, 2)))+c[1]*((diff(x[1](t), t))-(diff(x[0](t), t)))+k[1]*(x[1](t)-x[0](t)) = 0], [m[1]*(diff(diff(x[1](t), t), t))+c[1]*((diff(x[1](t), t)... (1.1)

m[1]は一階の質量です。モデリングを簡単にするためにベクトル形式にします。
x[1](t)は時間tの一階のx座標(横方向)です。
c[1]は粘性減衰定数です。k[1]は一階と地面の横のバネ定数です。
x[0](t)は地面のx軸方向座標で自身の横揺れの入力をこの関数で表現します。

同様に2階の運動方程式です。

> eq[2]:=m[2]*diff(x[2](t),t$2)+c[2]*(diff(x[2](t),t)-diff(x[1](t),t))+k[2]*(x[2](t)-x[1](t))=0;

(Typesetting:-mprintslash)([eq[2] := m[2]*(diff(x[2](t), `$`(t, 2)))+c[2]*((diff(x[2](t), t))-(diff(x[1](t), t)))+k[2]*(x[2](t)-x[1](t)) = 0], [m[2]*(diff(diff(x[2](t), t), t))+c[2]*((diff(x[2](t), t)... (1.2)

同じく3階の方程式です。

> eq[3]:=m[3]*diff(x[3](t),t$2)+c[3]*(diff(x[3](t),t)-diff(x[2](t),t))+k[3]*(x[3](t)-x[2](t))=0;

(Typesetting:-mprintslash)([eq[3] := m[3]*(diff(x[3](t), `$`(t, 2)))+c[3]*((diff(x[3](t), t))-(diff(x[2](t), t)))+k[3]*(x[3](t)-x[2](t)) = 0], [m[3]*(diff(diff(x[3](t), t), t))+c[3]*((diff(x[3](t), t)... (1.3)

この連立方程式は数式的に解いても非常に複雑になりますので意味がありません。
各定数にモデルの定数を具体的に設定し簡単化し数値的なシミュレーションを行います。

まず質量を決めます。
> m[1]:=100000;m[2]:=100000;m[3]:=100000;

(Typesetting:-mprintslash)([m[1] := 100000], [100000])
(Typesetting:-mprintslash)([m[2] := 100000], [100000])
(Typesetting:-mprintslash)([m[3] := 100000], [100000]) (1.4)

次に粘性減衰定数です。
> c[1]:=10;c[2]:=10;c[3]:=10;

(Typesetting:-mprintslash)([c[1] := 10], [10])
(Typesetting:-mprintslash)([c[2] := 10], [10])
(Typesetting:-mprintslash)([c[3] := 10], [10]) (1.5)

各階の間の横のバネ定数です。
> k[1]:=1000000000;k[2]:=1000000000;k[3]:=1000000000;

(Typesetting:-mprintslash)([k[1] := 1000000000], [1000000000])
(Typesetting:-mprintslash)([k[2] := 1000000000], [1000000000])
(Typesetting:-mprintslash)([k[3] := 1000000000], [1000000000]) (1.6)

地面の動きの入力の関数を定義します。
> x[0]:=t->sin(5*t)/10;

(Typesetting:-mprintslash)([x[0] := proc (t) options operator, arrow; 1/10*sin(5*t) end proc], [proc (t) options operator, arrow; 1/10*sin(5*t) end proc]) (1.7)

定数を決めると各階の運動方程式は次のように簡単になります。
> eq[1];

(Typesetting:-mprintslash)([100000*(diff(x[1](t), `$`(t, 2)))+10*(diff(x[1](t), t))-5*cos(5*t)+1000000000*x[1](t)-100000000*sin(5*t) = 0], [100000*(diff(diff(x[1](t), t), t))+10*(diff(x[1](t), t))-5*c...
(Typesetting:-mprintslash)([100000*(diff(x[1](t), `$`(t, 2)))+10*(diff(x[1](t), t))-5*cos(5*t)+1000000000*x[1](t)-100000000*sin(5*t) = 0], [100000*(diff(diff(x[1](t), t), t))+10*(diff(x[1](t), t))-5*c...
(1.8)

> eq[2];

(Typesetting:-mprintslash)([100000*(diff(x[2](t), `$`(t, 2)))+10*(diff(x[2](t), t))-10*(diff(x[1](t), t))+1000000000*x[2](t)-1000000000*x[1](t) = 0], [100000*(diff(diff(x[2](t), t), t))+10*(diff(x[2](...
(Typesetting:-mprintslash)([100000*(diff(x[2](t), `$`(t, 2)))+10*(diff(x[2](t), t))-10*(diff(x[1](t), t))+1000000000*x[2](t)-1000000000*x[1](t) = 0], [100000*(diff(diff(x[2](t), t), t))+10*(diff(x[2](...
(1.9)

> eq[3];

(Typesetting:-mprintslash)([100000*(diff(x[3](t), `$`(t, 2)))+10*(diff(x[3](t), t))-10*(diff(x[2](t), t))+1000000000*x[3](t)-1000000000*x[2](t) = 0], [100000*(diff(diff(x[3](t), t), t))+10*(diff(x[3](...
(Typesetting:-mprintslash)([100000*(diff(x[3](t), `$`(t, 2)))+10*(diff(x[3](t), t))-10*(diff(x[2](t), t))+1000000000*x[3](t)-1000000000*x[2](t) = 0], [100000*(diff(diff(x[3](t), t), t))+10*(diff(x[3](...
(1.10)

初期値問題として初期値を決めます。
> init:={x[1](0)=0,x[2](0)=0,x[3](0)=0};

(Typesetting:-mprintslash)([init := {x[2](0) = 0, x[1](0) = 0, x[3](0) = 0}], [{x[2](0) = 0, x[1](0) = 0, x[3](0) = 0}]) (1.11)

まず各階の初期(T=0)での位置です。
> init1:={D(x[1])(0)=0,D(x[2])(0)=0,D(x[3])(0)=0};

(Typesetting:-mprintslash)([init1 := {(D(x[1]))(0) = 0, (D(x[2]))(0) = 0, (D(x[3]))(0) = 0}], [{(D(x[1]))(0) = 0, (D(x[2]))(0) = 0, (D(x[3]))(0) = 0}]) (1.12)

Dは微分オペレータでこれは各階の初速すなわちt=0での速度を与えています。
これでdsolveにより解きます。unionは集合の和を行う演算で連率方程式と初期値を合わせています。

微分方程式の数値解を求めますがそのサンプリング点の時系列sampを生成します。

まず時間刻みを決めます。
> Delta[T] := 0.3e-2

(Typesetting:-mprintslash)([Delta[T] := 0.3e-2], [0.3e-2]) (1.13)

サンプリング点数(+1)mmを決めます。
> mm := 100

(Typesetting:-mprintslash)([mm := 100], [100]) (1.14)

サンプリング点の数列を生成します。
> samp := array([seq(Delta[T]*j, j = 0 .. mm)])

(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(Typesetting:-mprintslash)([samp := vector([0., 0.3e-2, 0.6e-2, 0.9e-2, 0.12e-1, 0.15e-1, 0.18e-1, 0.21e-1, 0.24e-1, 0.27e-1, 0.30e-1, 0.33e-1, 0.36e-1, 0.39e-1, 0.42e-1, 0.45e-1, 0.48e-1, 0.51e-1, 0....
(1.15)

求める関数の集合を設定します。
> funcs := {x[2](t), x[1](t), x[3](t)}

(Typesetting:-mprintslash)([funcs := {x[2](t), x[1](t), x[3](t)}], [{x[2](t), x[1](t), x[3](t)}]) (1.16)

dsolveで数値解を求めます。
> ans:=dsolve({eq[1],eq[2],eq[3]} union init union init1,funcs,numeric,output=samp):

解に保存されているデータは以下のようになっています。
> print(ans[1, 1])

(Typesetting:-mprintslash)([vector([t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t)])], [ans[1, 1]]) (1.17)

一階の位置x[1]の階の最初の部分を時間といっしょにプリントしてみます。
> for j from 0 to 10 do print(ans[2, 1][j+1, 1], ans[2, 1][j+1, 2]) end do; -1
for j from 0 to 10 do print(ans[2, 1][j+1, 1], ans[2, 1][j+1, 2]) end do; -1
for j from 0 to 10 do print(ans[2, 1][j+1, 1], ans[2, 1][j+1, 2]) end do; -1

0., 0.

0.300000000000000006e-2, 0.223987652748355613e-4

0.600000000000000014e-2, 0.176779860296743494e-3

0.899999999999999932e-2, 0.583305987424735944e-3

0.120000000000000002e-1, 0.133955571868308557e-2

0.149999999999999994e-1, 0.251177756100701187e-2

0.179999999999999988e-1, 0.412894198002203631e-2

0.210000000000000014e-1, 0.618012338469219742e-2

0.240000000000000004e-1, 0.861544601564112264e-2

0.269999999999999996e-1, 0.113505097173965657e-1

0.299999999999999990e-1, 0.142738983886376552e-1 (1.18)

プロットのパッケージを呼び出しておきます。
> with(plots); -1; with(plottools); -1

地面の動きと一階の動きをプロットしてみます。
> pp:=seq([ans[2,1][j+1,1],ans[2,1][j+1,2]],j=0..mm):
> display(plot(x[0](t),t=0..0.3),plots[pointplot]([pp],style=line));
Plot
次は二階の動きです。
> pp := seq([ans[2, 1][j+1, 1], ans[2, 1][j+1, 4]], j = 0 .. mm); -1
pp := seq([ans[2, 1][j+1, 1], ans[2, 1][j+1, 4]], j = 0 .. mm); -1
> display(plot(x[0](t), t = 0 .. .3), plots[pointplot]([pp], style = line))
display(plot(x[0](t), t = 0 .. .3), plots[pointplot]([pp], style = line))
Plot
次に三階の動きをプロットしてみます。
> pp := seq([ans[2, 1][j+1, 1], ans[2, 1][j+1, 6]], j = 0 .. mm); -1
pp := seq([ans[2, 1][j+1, 1], ans[2, 1][j+1, 6]], j = 0 .. mm); -1
> display(plot(x[0](t), t = 0 .. .3), plots[pointplot]([pp], style = line))
display(plot(x[0](t), t = 0 .. .3), plots[pointplot]([pp], style = line))
Plot
建物の動きをアニメーション表示してみましょう。

基準点の青のオブジェクトを定義します。
> cx:=cuboid([-0.1,1,0],[0.1,1.2,-0.1],color=blue):

地面のオブジェクトです。
> c0:=cuboid([-1,-1,0],[1,1,-0.1]):

各階のオブジェクトを用意します。
> c1:=cuboid([-0.5,-0.5,0.4],[0.5,0.5,0.5],color=gray):
> c2:=cuboid([-0.5,-0.5,0.9],[0.5,0.5,1.0],color=gray):
> c3:=cuboid([-0.5,-0.5,1.4],[0.5,0.5,1.5],color=gray):

まずはこれらを表示してみましょう。
> display({cx,c0,c1,c2,c3},scaling=constrained);
Plot
これらのオブジェクトを動かして建物の動きを表現します。

例によってフレームを保存する変数pを空にセットしておきます。
> p:=[];

(Typesetting:-mprintslash)([p := []], [[]]) (1.19)

時間を変えながらフレームを作成します。このとき解の関数に従って
オブジェクトをtranslateで動かします。
> for j from 0 to mm do
t:=ans[2,1][j+1,1]:
xx[0]:=x[0](t):
xx[1]:=ans[2,1][j+1,2]:
xx[2]:=ans[2,1][j+1,4]:
xx[3]:=ans[2,1][j+1,6]:
 pp:=cx:
 pp:=pp,translate(c0,xx[0](t),0,0):
 pp:=pp,translate(c1,xx[1](t),0,0):
 pp:=pp,translate(c2,xx[2](t),0,0):
 pp:=pp,translate(c3,xx[3](t),0,0):
 p:=p,display(pp):
end:
> display(p,insequence=true,scaling=constrained);
Plot

10階の建物の横揺れのシミュレーション

今度は10階の建物(10質点)の横揺れのモデルです。

またMapleを初期化します。
> restart:
プロットのパッケージを呼び出しておきます。
> with(plots):with(plottools):
forループを使いますので階数を表す変数nを10階にセットしておきます。
> n:=10;

(Typesetting:-mprintslash)([n := 10], [10]) (2.1)

十個(10階)の方程式をforループで一気に作成します。
> for j from 1 to n do
 eq[j]:=m[j]*diff(x[j](t),t$2)+c[j]*(diff(x[j](t),t)-diff(x[j-1](t),t))+k[j]*(x[j](t)-x[j-1](t))=0;
end;

(Typesetting:-mprintslash)([eq[1] := m[1]*(diff(x[1](t), `$`(t, 2)))+c[1]*((diff(x[1](t), t))-(diff(x[0](t), t)))+k[1]*(x[1](t)-x[0](t)) = 0], [m[1]*(diff(diff(x[1](t), t), t))+c[1]*((diff(x[1](t), t)...

(Typesetting:-mprintslash)([eq[2] := m[2]*(diff(x[2](t), `$`(t, 2)))+c[2]*((diff(x[2](t), t))-(diff(x[1](t), t)))+k[2]*(x[2](t)-x[1](t)) = 0], [m[2]*(diff(diff(x[2](t), t), t))+c[2]*((diff(x[2](t), t)...

(Typesetting:-mprintslash)([eq[3] := m[3]*(diff(x[3](t), `$`(t, 2)))+c[3]*((diff(x[3](t), t))-(diff(x[2](t), t)))+k[3]*(x[3](t)-x[2](t)) = 0], [m[3]*(diff(diff(x[3](t), t), t))+c[3]*((diff(x[3](t), t)...

(Typesetting:-mprintslash)([eq[4] := m[4]*(diff(x[4](t), `$`(t, 2)))+c[4]*((diff(x[4](t), t))-(diff(x[3](t), t)))+k[4]*(x[4](t)-x[3](t)) = 0], [m[4]*(diff(diff(x[4](t), t), t))+c[4]*((diff(x[4](t), t)...

(Typesetting:-mprintslash)([eq[5] := m[5]*(diff(x[5](t), `$`(t, 2)))+c[5]*((diff(x[5](t), t))-(diff(x[4](t), t)))+k[5]*(x[5](t)-x[4](t)) = 0], [m[5]*(diff(diff(x[5](t), t), t))+c[5]*((diff(x[5](t), t)...

(Typesetting:-mprintslash)([eq[6] := m[6]*(diff(x[6](t), `$`(t, 2)))+c[6]*((diff(x[6](t), t))-(diff(x[5](t), t)))+k[6]*(x[6](t)-x[5](t)) = 0], [m[6]*(diff(diff(x[6](t), t), t))+c[6]*((diff(x[6](t), t)...

(Typesetting:-mprintslash)([eq[7] := m[7]*(diff(x[7](t), `$`(t, 2)))+c[7]*((diff(x[7](t), t))-(diff(x[6](t), t)))+k[7]*(x[7](t)-x[6](t)) = 0], [m[7]*(diff(diff(x[7](t), t), t))+c[7]*((diff(x[7](t), t)...

(Typesetting:-mprintslash)([eq[8] := m[8]*(diff(x[8](t), `$`(t, 2)))+c[8]*((diff(x[8](t), t))-(diff(x[7](t), t)))+k[8]*(x[8](t)-x[7](t)) = 0], [m[8]*(diff(diff(x[8](t), t), t))+c[8]*((diff(x[8](t), t)...

(Typesetting:-mprintslash)([eq[9] := m[9]*(diff(x[9](t), `$`(t, 2)))+c[9]*((diff(x[9](t), t))-(diff(x[8](t), t)))+k[9]*(x[9](t)-x[8](t)) = 0], [m[9]*(diff(diff(x[9](t), t), t))+c[9]*((diff(x[9](t), t)...
 

(Typesetting:-mprintslash)([eq[10] := m[10]*(diff(x[10](t), `$`(t, 2)))+c[10]*((diff(x[10](t), t))-(diff(x[9](t), t)))+k[10]*(x[10](t)-x[9](t)) = 0], [m[10]*(diff(diff(x[10](t), t), t))+c[10]*((diff(x...
(Typesetting:-mprintslash)([eq[10] := m[10]*(diff(x[10](t), `$`(t, 2)))+c[10]*((diff(x[10](t), t))-(diff(x[9](t), t)))+k[10]*(x[10](t)-x[9](t)) = 0], [m[10]*(diff(diff(x[10](t), t), t))+c[10]*((diff(x...
(2.2)

また各定数を決めて方程式を簡単にします。
> for j to n do
 m[j]:=5000:
end:
粘性定数
> for j to n do
 c[j]:=1000:
end:
バネ定数
> for j to n do
 k[j]:=1.25*10^8:
end:
地面の動き(地震)。わずかに動いた状態(ステップ応答)。0がもとの位置。
> xx[0]:=t->0.00001;
(Typesetting:-mprintslash)([xx[0] := proc (t) options operator, arrow; 0.1e-4 end proc], [proc (t) options operator, arrow; 0.1e-4 end proc]) (2.3)
> x[0]:=xx[0];

(Typesetting:-mprintslash)([x[0] := xx[0]], [xx[0]]) (2.4)

ここでx[0]は関数となっています。
evalfで方程式を数値化します。
> for j to n do
 eq[j]:=evalf(eq[j]):
end;

(Typesetting:-mprintslash)([eq[1] := 5000.*(diff(x[1](t), `$`(t, 2)))+1000.*(diff(x[1](t), t))+125000000.0*x[1](t)-1250.000000 = 0.], [5000.*(diff(diff(x[1](t), t), t))+1000.*(diff(x[1](t), t))+125000...
(Typesetting:-mprintslash)([eq[1] := 5000.*(diff(x[1](t), `$`(t, 2)))+1000.*(diff(x[1](t), t))+125000000.0*x[1](t)-1250.000000 = 0.], [5000.*(diff(diff(x[1](t), t), t))+1000.*(diff(x[1](t), t))+125000...

(Typesetting:-mprintslash)([eq[2] := 5000.*(diff(x[2](t), `$`(t, 2)))+1000.*(diff(x[2](t), t))-1000.*(diff(x[1](t), t))+125000000.0*x[2](t)-125000000.0*x[1](t) = 0.], [5000.*(diff(diff(x[2](t), t), t)...
(Typesetting:-mprintslash)([eq[2] := 5000.*(diff(x[2](t), `$`(t, 2)))+1000.*(diff(x[2](t), t))-1000.*(diff(x[1](t), t))+125000000.0*x[2](t)-125000000.0*x[1](t) = 0.], [5000.*(diff(diff(x[2](t), t), t)...

(Typesetting:-mprintslash)([eq[3] := 5000.*(diff(x[3](t), `$`(t, 2)))+1000.*(diff(x[3](t), t))-1000.*(diff(x[2](t), t))+125000000.0*x[3](t)-125000000.0*x[2](t) = 0.], [5000.*(diff(diff(x[3](t), t), t)...
(Typesetting:-mprintslash)([eq[3] := 5000.*(diff(x[3](t), `$`(t, 2)))+1000.*(diff(x[3](t), t))-1000.*(diff(x[2](t), t))+125000000.0*x[3](t)-125000000.0*x[2](t) = 0.], [5000.*(diff(diff(x[3](t), t), t)...

(Typesetting:-mprintslash)([eq[4] := 5000.*(diff(x[4](t), `$`(t, 2)))+1000.*(diff(x[4](t), t))-1000.*(diff(x[3](t), t))+125000000.0*x[4](t)-125000000.0*x[3](t) = 0.], [5000.*(diff(diff(x[4](t), t), t)...
(Typesetting:-mprintslash)([eq[4] := 5000.*(diff(x[4](t), `$`(t, 2)))+1000.*(diff(x[4](t), t))-1000.*(diff(x[3](t), t))+125000000.0*x[4](t)-125000000.0*x[3](t) = 0.], [5000.*(diff(diff(x[4](t), t), t)...

(Typesetting:-mprintslash)([eq[5] := 5000.*(diff(x[5](t), `$`(t, 2)))+1000.*(diff(x[5](t), t))-1000.*(diff(x[4](t), t))+125000000.0*x[5](t)-125000000.0*x[4](t) = 0.], [5000.*(diff(diff(x[5](t), t), t)...
(Typesetting:-mprintslash)([eq[5] := 5000.*(diff(x[5](t), `$`(t, 2)))+1000.*(diff(x[5](t), t))-1000.*(diff(x[4](t), t))+125000000.0*x[5](t)-125000000.0*x[4](t) = 0.], [5000.*(diff(diff(x[5](t), t), t)...

(Typesetting:-mprintslash)([eq[6] := 5000.*(diff(x[6](t), `$`(t, 2)))+1000.*(diff(x[6](t), t))-1000.*(diff(x[5](t), t))+125000000.0*x[6](t)-125000000.0*x[5](t) = 0.], [5000.*(diff(diff(x[6](t), t), t)...
(Typesetting:-mprintslash)([eq[6] := 5000.*(diff(x[6](t), `$`(t, 2)))+1000.*(diff(x[6](t), t))-1000.*(diff(x[5](t), t))+125000000.0*x[6](t)-125000000.0*x[5](t) = 0.], [5000.*(diff(diff(x[6](t), t), t)...

(Typesetting:-mprintslash)([eq[7] := 5000.*(diff(x[7](t), `$`(t, 2)))+1000.*(diff(x[7](t), t))-1000.*(diff(x[6](t), t))+125000000.0*x[7](t)-125000000.0*x[6](t) = 0.], [5000.*(diff(diff(x[7](t), t), t)...
(Typesetting:-mprintslash)([eq[7] := 5000.*(diff(x[7](t), `$`(t, 2)))+1000.*(diff(x[7](t), t))-1000.*(diff(x[6](t), t))+125000000.0*x[7](t)-125000000.0*x[6](t) = 0.], [5000.*(diff(diff(x[7](t), t), t)...

(Typesetting:-mprintslash)([eq[8] := 5000.*(diff(x[8](t), `$`(t, 2)))+1000.*(diff(x[8](t), t))-1000.*(diff(x[7](t), t))+125000000.0*x[8](t)-125000000.0*x[7](t) = 0.], [5000.*(diff(diff(x[8](t), t), t)...
(Typesetting:-mprintslash)([eq[8] := 5000.*(diff(x[8](t), `$`(t, 2)))+1000.*(diff(x[8](t), t))-1000.*(diff(x[7](t), t))+125000000.0*x[8](t)-125000000.0*x[7](t) = 0.], [5000.*(diff(diff(x[8](t), t), t)...

(Typesetting:-mprintslash)([eq[9] := 5000.*(diff(x[9](t), `$`(t, 2)))+1000.*(diff(x[9](t), t))-1000.*(diff(x[8](t), t))+125000000.0*x[9](t)-125000000.0*x[8](t) = 0.], [5000.*(diff(diff(x[9](t), t), t)...
(Typesetting:-mprintslash)([eq[9] := 5000.*(diff(x[9](t), `$`(t, 2)))+1000.*(diff(x[9](t), t))-1000.*(diff(x[8](t), t))+125000000.0*x[9](t)-125000000.0*x[8](t) = 0.], [5000.*(diff(diff(x[9](t), t), t)...

(Typesetting:-mprintslash)([eq[10] := 5000.*(diff(x[10](t), `$`(t, 2)))+1000.*(diff(x[10](t), t))-1000.*(diff(x[9](t), t))+125000000.0*x[10](t)-125000000.0*x[9](t) = 0.], [5000.*(diff(diff(x[10](t), t...
(Typesetting:-mprintslash)([eq[10] := 5000.*(diff(x[10](t), `$`(t, 2)))+1000.*(diff(x[10](t), t))-1000.*(diff(x[9](t), t))+125000000.0*x[10](t)-125000000.0*x[9](t) = 0.], [5000.*(diff(diff(x[10](t), t...
(2.5)

初期値の準備。各階の初期位置を集合で与えます。全て0とします。
> init:={};

(Typesetting:-mprintslash)([init := {}], [{}]) (2.6)

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

各階の初速を与えます。全て0とします。
> init1:={};

(Typesetting:-mprintslash)([init1 := {}], [{}]) (2.7)

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

連立方程式の集合をセットします。
> equ:={};

(Typesetting:-mprintslash)([equ := {}], [{}]) (2.8)

> for j to n do
 equ:=equ union {eq[j]}
end:

解く解の関数の集合をセット。
> funcs:={};

(Typesetting:-mprintslash)([funcs := {}], [{}]) (2.9)

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

結果を保存するサンプリング時間とサンプル数の設定をします。
> d:=0.002;l:=50;

(Typesetting:-mprintslash)([d := 0.2e-2], [0.2e-2])
(Typesetting:-mprintslash)([l := 50], [50]) (2.10)

0.001秒毎。100サンプル(100回)計算

数値解の結果を求めるサンプル時間を生成します。繰り返しにseqコマンドを利用します。
> samp:=array([seq(d*j,j=0..l)]):
10個の連立の運動方程式をdsolveで解きます。このとき数値解を求めますのでnumericを指定します。
outputでサンプル時間を指定しサンプル点での解を求めます。
> ans:=dsolve(equ union init union init1,funcs,numeric,output=samp):
output=で指定したサンプル時間にたいしての数値解が得られます。
得られるデータのリストです。各階の位置と速度(一階微分)。繰り返しにはseqコマンドが便利です。
> seq(ans[1,1][j],j=1..2*n+1);

(Typesetting:-mprintslash)([t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t)...
(Typesetting:-mprintslash)([t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t)...
(Typesetting:-mprintslash)([t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t)...
(2.11)

一階の位置をプロットします。サンプルデータ(リスト構造)はlistplotでプロットできます。
> listplot([seq(ans[2,1][j,2],j=1..l)]);
Plot
二階の位置をプロット。
> listplot([seq(ans[2,1][j,4],j=1..l)]);
Plot
10階の位置をプロットします。
> listplot([seq(ans[2,1][j,20],j=1..l)]);
Plot
またアニメーションを作成しましょう。
まず各階を表すオブジェクトを作成しましょう。
> cc:=cuboid([-1,-1,0],[1,1,0.2]):
> display(cc,scaling=constrained,color=gray);
Plot
フレームを保存する変数を準備。空にしておきます。
> p:=[];

(Typesetting:-mprintslash)([p := []], [[]]) (2.12)

時間毎の各階の位置をフレームに書き込みます。
動きはわかりやすいように実際の値よりは拡大(400倍)して動かしています。
> for j from l/2  to 7*l/8 do
 t:=ans[2,1][j+1,1]:
 pp:=cuboid([-2,-2,0],[2,2,-0.1]):
 for i to n do
   pp:=pp,translate(cc,400*ans[2,1][j+1,i*2],0,i):
 end:
 p:=p,display(pp):
end:

生成されたフレームをアニメーション表示します。

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